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AOUADA,  DJAMILA.  Geometric,  Statistical,  and  Topological  Modeling  of  Intrinsic  Data 
Manifolds:  Application  to  3D  Shapes.  (Under  the  direction  of  Professor  Hamid  Krim). 

The  increasing  size  and  complexity  of  data  often  invokes  the  extraction  of  infor¬ 
mation  from  their  reduced  representations  while  preserving  their  inherent  structure.  In  this 
thesis,  we  explore  the  statistical,  geometric  and  topological  intrinsic  information  contained 
in  high  dimensional  data.  We  focus  on  applications  related  to  3-dimensional  objects,  and 
model  their  2-dimensional  surfaces  using  compact  curved-skeletal  models  that  we  refer  to 
as  “squigraphs” .  These  models  are  multi-level  representations  that  superpose  global  topo¬ 
logical  and  local  geometric  3D  shape  descriptors.  Squigraphs  are  subsequently  used  for 
classification,  and  ensure  a  high  discrimination  between  in-class  3-dimensional  shapes. 

The  extraction  of  squigraphs  starts  by  sampling  the  surface  of  an  object  for  a 
resulting  set  of  curves.  This  may  be  accomplished  by  defining  an  appropriate  intrinsic 
characteristic  function  on  the  surface  itself,  referred  to  as  a  Morse  function;  which  we  use 
in  a  two-phase  approach.  To  ensure  the  invariance  of  the  final  representation  to  isometric 
transforms,  we  choose  the  Morse  function  to  be  an  intrinsic  global  geodesic  function.  The 
first  phase  is  a  coarse  representation  through  a  reduced  topological  Reeb  graph.  We  use  it  for 
a  meaningful  decomposition  of  shapes  into  primitives.  At  the  second  phase,  we  add  detailed 
geometric  information  by  tracking  the  evolution  of  Morse  function’s  level  curves  along  each 
primitive.  We  then  embed  the  manifold  corresponding  to  this  evolution  of  curves  into  R3, 
and  obtain  a  simple  space  curve.  We  further  define  a  Riemannian  metric  to  quantitatively 
compare  the  geometry  of  shapes. 

We  point  the  flexibility  of  our  techniques  for  other  applications,  namely,  face  recog¬ 
nition,  behavioral  modeling,  and  sensor  network  data  analysis.  While  all  these  applications 
face  the  same  curse  of  dimensionality,  we  show  that  they  may  be  formalized  under  similar 
geometrical  settings. 
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Chapter  1 

Introduction 


The  last  decade  has  seen  an  important  growth  in  applied  problems  that  are  tightly 
related  to  the  explosive  increase  in  data  sizes  and  complexities.  The  nature  of  these  appli¬ 
cations  calls  for  a  parsimonious  representation  of  the  data.  The  most  challenging  constraint 
associated  with  this  task  is  the  preservation  of  the  inherent  data  structure  in  the  course  of 
analysis  and  modeling. 

In  this  thesis,  we  explore  the  statistical,  geometric  and  topological  intrinsic  information 
contained  in  high  dimensional  data  structures.  Our  objective  is  to  define  simpler  models 
with  significantly  reduced  complexities,  which  will  subsequently  facilitate  the  usual  compu¬ 
tationally  demanding  applications  such  as  data  classification  and  data  retrieval.  Reducing 
the  dimension  of  a  data  set  is  often  based  on  information  that  is  characteristic  and  restricted 
to  a  small  portion  of  the  whole  data  space  to  be  fully  expressed  or  summarized.  Finding 
and  correctly  manipulating  the  necessary  and  sufficient  features  constitute  the  central  goal 
of  our  work. 

1.1  Motivations  and  overview 

In  the  present  thesis,  we  focus  on  modeling  the  2-dimensional  surface  of  a  3D  ob¬ 
ject.  Our  interest  in  this  kind  of  data  is  motivated  by  a  variety  of  applications  ranging  from 
target  recognition  and  face  identification  to  shape  retrieval  and  behavior  analysis. 

While  a  number  of  available  3D  shape  modeling  techniques  [1,  2,  3,  4,  5,  6]  yield  satisfactory 
object  classification  results,  many  applications  require  a  more  refined  and  efficient  identi- 
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fication/recognition  of  objects  in  a  given  class.  To  that  end,  we  invoke  Morse  theory  in  a 
two-phase  approach.  To  ensure  the  invariance  of  a  3D  object  final  representation  to  isomet¬ 
ric  transforms,  we  choose  the  Morse  function  to  be  a  simple  and  intrinsic  global  geodesic 
function  defined  on  its  surface.  At  first,  a  coarse  representation  through  a  reduced  topolog¬ 
ical  Reeb  graph  is  obtained.  This  is  in  turn  used  in  a  meaningful  decomposition  of  shapes 
into  primitives.  Following  this  step,  we  add  detailed  geometric  information  by  tracking  the 
evolution  of  Morse  function’s  level  curves  along  each  primitive.  We  subsequently  embed  the 
manifold  corresponding  to  the  evolution  of  these  curves  into  M3,  and  obtain  a  simple  space 
curve.  By  combining  phases  one  and  two,  we  construct  new  graphs  rich  in  topological  and 
geometric  information  that  we  refer  to  as  squigraphs. 

Our  proposed  analysis  framework  essentially  exploits  curves  in  space  (which  are 
representative  of  3D  objects).  As  we  further  elaborate  later,  spatial  curves  may  indeed 
be  exploited  in  various  configurations.  They  may,  for  instance,  be  extracted  as  contours 
of  landmark  surfaces  [7],  as  level  curves  of  a  Morse  function  [8],  or  also  as  elements  of 
curved  skeletons  [9,  10].  All  these  techniques,  despite  their  differences,  commonly  rely 
on  curves’  properties  in  solving  computer  vision  problems.  This  approach  is  conceptually 
motivated  by  the  fact  that  curves  in  3D  are  fairly  well  known  geometric  entities;  moreover, 
under  some  conditions,  they  can  accurately  describe  the  overall  geometry  of  an  object  in 
3D  space  [7].  Translating  the  constraints  of  3D  shape  representation  techniques  to  curves, 
reduces  the  level  of  analytical  difficulty  associated  with  the  3D  representation  problem  and 
makes  it  more  tractable.  To  this  end,  we  focus  our  efforts  on  defining  an  effective  modeling 
framework  for  curves  in  space  along  with  the  corresponding  Riemannian  metric. 

While  3D  objects  form  the  core  of  the  applications  of  interest  in  this  thesis,  we  also 
discuss  the  exploration  of  other  types  of  complex  data  ( e.g video  frames,  sensor  network 
data).  Our  techniques  are  thus  readily  adaptable  to  various  applications  as  long  as  the 
measured  data  may  be  correctly  assumed  to  lie  on  a  smooth  manifold.  Indeed,  we  show 
how  our  proposed  Whitney  modeling  of  curves  can  be  used  to  represent  the  evolution  of 
video  frames,  where  a  frame  is  equated  to  a  segmented  silhouette.  The  evolution  of  a  set  of 
silhouettes  that  is  ordered  in  time  becomes  a  1-dimensional  manifold  that  can  be  embedded 
in  3D. 

A  less  visually  intuitive  problem  is  that  of  analyzing  sensor  network  data.  Our  goal 
is  to  develop  mathematical  models  which  will  lead  to  a  clear  understanding  of  the  interde- 


3 


pendencies  between  different  nodes  within  one  network,  as  well  as  of  the  interdependencies 
between  different  networks.  Such  models  will  be  crucial  in  guiding  the  implementation  of 
network  infrastructures/protocols  which  are  robust  and  resilient  to  intrusions.  Our  approach 
starts  by  abstracting  networks  from  their  physical  layered  description.  We  view  a  network 
as  a  set  of  features  which  best  characterize  it.  The  same  feature  space  may  house  multi¬ 
ple  distinct  networks.  Our  hypothesis  is  that  each  network  is  characterized  by  a  smooth 
topologically  homogeneous  manifold  embedded  in  the  network  feature  space;  therefore,  we 
propose  to  carry  out  our  work  by  relying  on  notions  from  topological  and  geometric  mani¬ 
fold  learning.  The  interdependencies  between  networks  may  thus  be  reflected  by  statistical 
correlations  between  the  different  space  elements  (e.g.,  manifolds  and  points),  where  data 
points  on  a  manifold  are  representative  of  nodes  in  a  network.  In  addition,  an  attack  or  any 
intrusion  is  interpreted  as  having  the  effect  of  destroying  the  initial  topological  structure 
of  a  manifold.  Our  work  is  hence  geared  toward  detecting  and  determining  any  topological 
changes. 

1.2  Contributions 

We  summarize  our  main  contributions  bellow. 

•  We  have  investigated  the  representation,  classification  and  recognition  of  3D  objects. 
We  defined  and  implemented  a  complete  machinery  for  a  precise  identification  and 
recognition  of  3D  objects  from  multiple  classes  and  within  one  class.  In  light  of  the 
fact  that  shapes  can  be  very  complex,  all  our  proposed  techniques  promote  the  neces¬ 
sity  to  distinguish  between  global  and  local  features,  or  between  the  coarse  and  fine 
representations.  For  each  of  our  proposed  applications,  we  first  determine  the  group 
of  transforms  to  which  the  corresponding  representations  are  to  be  invariant.  Our 
focus  has  been  on  isometries,  as  they  gather  rigid  and  non-rigid  shapes  and  cover  a 
wide  range  of  applications  from  biomedical  imaging  to  computer  animations.  To  that 
end,  we  make  an  extensive  use  of  Morse  theory  throughout  this  work.  Specifically, 
we  choose  to  define  a  fully  intrinsic  global  geodesic  Morse  function  on  the  surface 
of  3D  objects.  We  exploit  the  information  contained  in  this  function  in  a  two-phase 
approach.  The  first  phase  is  a  coarse  representation  through  a  reduced  topological 
Reeb  graph,  which  is  instrumental  in  decomposing  a  complex  shape  into  primitives. 
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In  a  second  phase,  the  objective  is  to  provide  a  more  fine  representation.  Taking 
advantage  of  the  same  Morse  function,  we  proceed  to  analyze  the  evolution  of  its  level 
curves  along  each  primitive  and  further  reduce  its  representation  to  a  simple  curve 
in  3D  by  using  the  Whitney  embedding  theorem.  In  a  very  ludic  way,  we  propose  to 
combine  phases  one  and  two  to  construct  “squigraphs”  [10].  Experiments  show  that 
squigraphs  are  more  general  than  existing  techniques  [11];  they  achieve  very  competi¬ 
tive  classification  rates  in  comparison  to  those  achieved  by  classical  shape  descriptors. 
Their  performance,  however,  becomes  clearly  superior  when  finer  classification  and 
identification  operations  are  desired. 

•  In  the  same  spirit  of  intrinsic  modeling,  we  introduce  a  Riemannian  metric  to  compare 
curves  in  space.  The  classification/recognition  of  3D  objects  by  way  of  squigraphs, 
forms  the  primary  motivation  of  developing  a  quantitative  comparison  framework  for 
space  curves.  Our  additional  objective  is  to  limit  the  invariance  of  curves’  shapes 
to  rigid  transforms.  Inspired  by  results  in  visual  cognition  and  human  perception, 
we  propose  a  technique  that  simultaneously  verifies  all  the  following  properties:  high 
discrimination  level,  independence  of  parametrization  and  independence  of  any  refer¬ 
ence  point,  uniqueness  property,  as  well  as  a  good  preservation  of  the  correspondence 
between  curves  to  effectively  tackle  partial  matching  problems  [12,  13]. 

•  In  a  different  context  yet  with  the  same  view  of  geometry-driven  analysis,  we  seek 
to  understand  network  data  through  statistical  manifold  learning.  The  key  idea  in 
this  problem  is  to  consider  data  samples  as  lying  on  one  or  more  manifolds.  From 
there,  one  may  start  studying  the  topo-geometric  and  statistical  properties  of  the  data 
by  calling  upon  well-established  theorems.  The  main  contribution  of  this  work  is  to 
harness  space  and  time  information  with  the  correlation  between  networks’  nodes. 

1.3  Outline 

The  present  thesis  is  organized  as  follows: 

•  We  cover  in  Chapter  2  the  necessary  background  to  understand  the  developments 
in  subsequent  chapters.  We  start  by  describing  the  common  discrete  representation 
of  3D  objects,  i.e.,  the  triangulated  meshes;  and  by  providing  a  brief  introduction 
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to  Morse  theory  along  with  a  Global  Geodesic  Function  (GGF)  as  a  Morse  function 
which  we  repeatedly  use  throughout  this  thesis. 

•  Chapter  3  presents  a  statistical  analysis  of  the  GGF  and  the  new  resulting  object 
classification  algorithm.  We  refer  to  this  algorithm  as  Classification  using  object 
Characteristic  Resolution  (CCR). 

•  Chapter  4  describes  the  partitioning  of  3D  objects  using  their  Reeb  graphs.  The 
extraction  of  Reeb  graphs  is  herein  based  on  level  curves  of  a  Morse  function  instead 
of  its  intervals. 

•  Chapter  5  defines  squigraphs  and  the  Whitney  embedding  for  the  geometric  modeling 
of  shapes. 

•  In  Chapter  6,  we  define  a  new  similarity  invariant  signature  that  provides  a  Rieman- 
nian  metric  for  an  intrinsic  comparison  of  shapes  via  their  squigraphs. 

•  Chapter  7  demonstrates,  through  extensive  experiments,  the  effectiveness  of  squigraphs 
when  combined  with  the  distance  metric  defined  in  Chapter  6. 

•  Chapter  8  illustrates  the  problem  of  sensor  network  analysis  with  the  same  spirit  of 
intrinsic  modeling. 

•  We  provide  conclusions  and  possible  future  research  directions  in  Chapter  9. 
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Chapter  2 

Background 


Research  interest  in  3D  object  analysis  has  witnessed  an  explosive  growth  over 
the  last  few  years.  While  this  may  in  part  be  explained  by  an  equally  impressive  growth 
in  computing  power,  the  availability  of  3D  data  acquisition  systems  and  ready  access  to  it 
at  relatively  low  cost,  have  been  key  in  addressing  the  numerous  problems  which  arise  in 
applications.  Laser  scanners,  ranging  cameras  and  others  have  indeed  made  solutions  to 
multimedia  applications,  biometrics,  and  computer  graphics  more  realistic  and  affordable. 
The  wide  distribution  of  3D  data  over  the  internet  at  no-cost  is  also  testimony  to  the  high 
level  of  interest  in  the  area  [14,  15,  16,  17,  18]. 

In  what  follows,  we  describe  the  discrete  representation  that  is  often  used  to  present  3D  ob¬ 
jects.  We  also  describe  the  continuous  surface  of  these  objects  that  we  use  in  our  theoretical 
formulations. 

2.1  Digital  and  continuous  3D  objects 

3D  objects  are  commonly  represented  as  polygonal  or  triangulated  meshes.  Two 
matrices  T  and  V  are  the  usual  digital  representation  of  a  3D  object,  and  are  often  of 
prohibitive  size.  This  in  turn,  unveils  the  computational  challenge  often  encountered  in 
practice  and  particularly  in  3D  shape  classification  and  recognition  applications.  A  tri¬ 
angulated  mesh  M  is  thus  denoted  by  M  —  (V,JE),  where  V  =  [vi,...,vn]  is  the  matrix 
of  vertices,  such  that  contains  the  Euclidean  coordinates  of  the  vertex  labeled  i,  i.e., 
\~i  —  [ Xi,yi,Zi]T  and  T  —  [fi,  ...,fm]  is  the  matrix  of  faces,  such  that  f )  contains  the  labels 
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Figure  2.1:  Digital  and  continuous  representations  of  a  3D  object. 


of  all  the  vertices  of  the  face  under  the  label  j,  j  =  1,  •  •  •  ,  m.  We  refer  to  the  number  n 
as  the  resolution  for  a  digital  representation  of  a  3D  object.  Mathematically,  however,  we 
view  3D  objets  as  smooth  and  compact  surfaces  S  embedded  in  M3.  We  show  in  Figure  2.1 
an  example  of  a  triangulated  mesh,  and  the  corresponding  continuous  surface. 

2.2  Morse  theory 

In  order  to  efficiently  represent  3D  objects  and  extract  their  topological  and  ge¬ 
ometric  features,  we  make  an  extensive  use  of  Morse  theory  [19,  37].  Morse  theory  states 
that  it  is  possible  to  define  a  particular  smooth  function  /(•)  on  a  smooth  surface  S  (or 
more  generally  a  manifold  <S),  and  track  its  critical  points  in  order  to  study  the  topology  of 
S.  Such  a  function  /(•)  is  called  a  Morse  function  and  is  defined  as  follows: 

Definition  1  (Morse  function)  A  smooth  function  f  :  S  — >  R  on  a  smooth  manifold  S 
is  called  Morse  if  all  of  its  critical  points  are  non- degenerate. 

A  standard  Morse  function  is  a  height  function  [37,  21,  22].  The  image  of  a  point  p  on 
S  via  the  height  function  is  reduced  to  its  z  coordinate  as  shown  in  Figure  2.2  (b).  Any 
representation  that  is  based  on  the  height  function  is  clearly  varying,  as  critical  points  of  a 


Figure  2.2:  Illustration  of  different  Morse  functions:  (a)  Initial  object,  (b)  height  function 
on  the  surface  of  a  double  torus,  (c)  GGF  function  on  the  surface.  (Best  visualized  in  color) 


surface  in  R3  will  change  as  a  result  of  its  mere  rotation.  Since  shape  descriptors  are,  for 
a  large  number  of  applications,  required  to  be  invariant  to  similarity  transforms1,  a  choice 
of  an  appropriate  Morse  function  is  critical.  In  [1],  Hilaga  et  al.  define  a  Morse  function 
for  generic  metrics  on  surfaces,  that  is  invariant  to  isometric  transforms.  Specifically,  this 
function  is  defined  at  every  point  v  on  S  as  the  integral  of  the  geodesic  distance  d(v,  p) 
from  v  to  all  other  points  p  on  5, 

/(V)  =  [  d(v,p)d<S.  (2.1) 

J  pe<S 

2.3  Global  Geodesic  Function  GGF 

Hilaga  et  aV s  discrete  approximation  of  Eq.  (2.1),  that  we  herein  refer  to  as  the 
first  approximation  /appr(*),  is  defined  as  follows: 

fappr ( v )  =  ^  d(v,  bj)  •  area(bj),  (2.2) 

i 

where  {b^=o,i,...  is  a  finitely  countable  set  of  base  vertices  scattered  on  S  and  area(b^) 
is  the  area  that  b^  occupies,  such  that,  JR  area(b i)  ~  area(5).  Also,  for  an  accurate 
result,  Hilaga  et  al  emphasize  the  need  for  a  mesh  preparation  through  two  operations: 
generation  of  short-cut  edges  based  on  a  manually  chosen  threshold  and  a  subdivision  of 
the  mesh.  Let’s  consider  the  object  “droplet”  of  Figure  2.3.  and  compute  its  integrated 


1  Similarity  transforms  include:  translation,  rotation,  and  scaling,  or  any  combination  thereof. 
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Figure  2.3:  Comparison  of  the  two  approximations  of  the  integrated  geodesic  Morse  func¬ 
tion.  (a)  and  (b)  are,  respectively,  the  front  and  bottom  views  of  the  distribution  of  fappr(m ) 
on  the  object  “droplet”,  (c)  and  (d)  are,  respectively,  the  front  and  bottom  views  of  the 
distribution  of  g(-)  on  the  object  “droplet”.  (Best  visualized  in  color) 


geodesic  function.  Because  of  the  perfect  symmetry  of  this  3D  shape  with  respect  to  the 
z  axis,  one  should  expect  the  function  /(•)  defined  in  Eq.  (2.1)  to  have  all  its  level  sets 
exactly  parallel  to  the  xy  plane,  i.e.,  horizontally  constant.  We  find,  however,  that  failure 
in  achieving  an  appropriate  prepossessing  of  the  mesh,  may  drastically  affect  the  distribution 
of  fappr ( ' ) ’  fhe  firsf  approximation  of  /(•)  as  defined  in  Eq.  (2.2).  In  Figure  2.3,  we  see  that 
the  colors  on  (a)  are  not  uniform  as  they  should  be  and  are  on  (c),  hence  our  motivation 
in  adopting  a  different  approximation  of  /(•)  as  our  Morse  function,  all  the  while  keeping 
similar  attractive  properties;  This  includes  its  full  invariance  to  isometric  transformations, 
and  a  further  improvement  of  the  robustness  to  surface  meshing  and  noise.  We  hence  use 
the  approximation  #(•),  that  we  first  defined  in  [23],  and  referred  to  as  the  Global  Geodesic 
Function  (GGF),  such  that: 


g(y) 


Epes  d(v>  P) 


maxqe<5 


(Xpe5d(<P  P)) 


(2.3) 


Unless  specified  otherwise,  this  is  the  GGF  used  throughout  this  thesis,  and  illustrated  in 
Figure  2.4,  and  in  Figure  2.2  (c),  where  the  surface  S  is  the  previously  seen  double  torus 
illustrated  in  Figure  2.2  (a).  In  addition  to  its  independence  of  any  reference  point,  an 
important  property  of  this  Morse  function  is  its  invariance  to  isometric  transformations 
(because  it  is  intrinsic),  thus  yielding  a  consistent  and  unique  characterization  of  an  object 
surface  as  illustrated  in  Figure  2.5.  The  independence  of  reference  points  is  achieved  by  the 
geodesic  integration  procedure  at  each  point.  The  normalization  of  the  functional  shown  in 
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Figure  2.4:  Examples  of  the  global 

geodesic  function  on  different  3D  objects. 


Figure  2.5:  Invariance  to  isometric  trans¬ 
formations:  (A)  Original  pose  and  a  cor¬ 
responding  PDF  of  its  GGF  (B)  Same 
subject  after  deformation  and  noise  ad¬ 
dition.  At  the  bottom  the  PDF  of  its 
GGF. 


Eq.  (2.3)  ensures  an  invariance  to  scaling  by  making  the  range  of  g(-)  coincide  with  [go,  1], 
with  go  —  max" ^5 ^P€<S  d(v  p)  an<^  In  Figure  2.2  (c),  we  show  the  color  coding  map  for 

g(-)  used  on  all  our  objects.  In  practice,  the  GGF  for  each  vertex  is  obtained  by  computing 
the  latter’s  geodesic  distance  to  all  other  vertices.  This  normalization  also  obviates  the 
explicit  computation  of  the  surface  element.  This  is  efficiently  realized  with  the  well  known 
Dijkstra  algorithm  whose  complexity  is  0(N2  log  N)  [24]. 
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Chapter  3 

Statistical  Analysis 


In  this  Chapter,  we  present  an  economical  and  fast  classification  algorithm  for  3D 
objects.  Classifying  3D  shapes  independently  of  any  rotation,  translation,  scaling  transfor¬ 
mation  and/or  non-elastic  deformation  is  of  great  importance  in  many  applications  includ¬ 
ing  multimedia,  communications,  computer  graphics  and  biometrics.  While  many  successful 
methods  exist,  they  are  of  relatively  high  complexity  for  the  task  of  putting  objects  into 
groups.  We  propose  a  classification  strategy  that  is  far  less  complex  and  still  very  robust 
to  pose  and  articulation  changes.  Our  method  is  based  on  the  GGF  defined  in  Section  2.3. 
Each  object  class  is  defined  by  only  four  parameters  obtained  during  a  learning  stage.  We 
use  these  parameters  in  the  decision  of  a  class  membership.  Experimental  results  demon¬ 
strate  low  cost,  efficiency,  and  robustness  to  resolution  and  data  benchmarks  of  the  proposed 
approach. 

3.1  Introduction 

We  propose  to  define  a  3D  object  classification  technique  where  we  characterize 
each  class  with  four  parameters  obtained  from  the  global  geodesic  shape  function.  The 
first  three  parameters  consist  of  a  variability  measure  and  two  characteristic  resolutions. 
The  characteristic  resolutions  indicate  how  to  maximally  and  effectively  reduce  the  overall 
computational  complexity  with  little  or  no  loss  of  information  about  a  given  object.  These 
parameters  are  generally  class  specific  and  may,  hence,  be  used  as  a  first  order  classifier.  A 
refined  level  of  discrimination  is  achieved  by  introducing  a  threshold  measure  as  a  fourth 
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Figure  3.1:  Illustration  of  the  visual  effect  of  the  resolution  reduction  on  the  GGF  of  a  3D 
object.  Resolution  is  decreasing  from  case  (a)  to  (f).  An  abrupt  change  in  the  GGF  occurs 
at  (d). 


parameter.  To  further  improve  on  the  discrimination  power  of  these  parameters,  we  propose 
a  subsequent  zooming  comparison  procedure  to  focus  on  distinctive  features  among  objects, 
and,  hence,  provide  a  potential  for  recognition1.  The  remainder  of  this  chapter  is  organized 
as  follows.  The  next  section  gives  the  formulation  of  the  problem  along  with  a  presentation 
of  the  key  tools  that  are  used.  The  detailed  algorithm  is  explained  in  Section  3.3.  The  sub¬ 
steps  are  illustrated  as  well.  In  Section  3.4,  experimental  results  that  aimed  to  demonstrate 
the  robustness  and  the  efficiency  of  the  method  are  given.  Finally,  Section  3.5  summarizes 
the  present  Chapter,  and  proposes  future  projects  as  a  continuity  to  the  present  statistical 
analysis. 

3.2  Problem  formulation 

Given  a  set  of  N  classes  of  objects  {Cj ,  C%, . . . ,  CV},  our  goal  is  to  decide  on  a 
class  membership  of  an  object  O.  To  efficiently  carry  through  such  a  task,  we  define  a  set 
of  feature  parameters  for  an  object  and  to  construct  a  systematic  procedure  for  2D  surface 
comparison. 

3.2.1  Resolution  evaluation 

As  stated  earlier,  the  computational  complexity  is  a  predominant  factor  in  3D 
object  processing.  We  have  experimentally  also  established  that  a  uniform  subdivision  of 


1  Recognition  here  means  a  more  refined  level  of  classification 
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Resolution  Evalution  with  the  highest  resolution  RQ  as  reference 
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Resolution  Evaluation  with  the  characteristic  resolution  R  as  reference 
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Figure  3.2:  Resolution  evaluation  for  a  dog  from  Princeton’s  benchmark. 


a  mesh  preserves  a  near  invariance  of  a  histogram  of  the  GGF,  for  different  resolutions 
of  an  object  relative  to  that  of  its  finest  resolution.  A  sharp  and  sudden  change  in  the 
overall  distribution  of  the  GGF  as  illustrated  in  Figure  3.1,  takes  place  and  the  visual 
perception  is  perturbed  when  the  resolution  is  too  low  to  exhibit  the  desired  invariance  of 
the  distribution.  The  reference  distribution  of  any  3D  object  in  a  class  of  interest  is  that 
of  its  GGF  computed  at  the  finest  resolution.  This  distribution  is  subsequently  compared 
to  those  obtained  by  progressively  decreasing  resolutions.  To  this  end,  we  use  the  Jensen- 
Shannon  Divergence  ( JSD )  as  a  distance  function  between  two  distributions  [3,  25].  For 
K  different  resolutions,  Rq  >  R\  >  . . .  >  Rk- i,  of  an  object  and  their  K  corresponding 
PDFs  P#. ,  i  =  0, 1, . . . ,  K  —  1,  the  JSD  distance  is  defined  as  follows 

JSD(PRl,PR2 )  =H(J2  \pRi)  -  E 

1= 1,2  Z  1= 1,2  Z 


(3-1) 
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where  H  is  the  Jensen- Shannon  entropy  defined  by: 

l=L 

H{P %)  =  -  PRi(l),  (3-2) 

1  =  1 

with  L  <  Ri  and  Pr.  (Z),  l  =  1,  2, . . . ,  L,  being  the  elements  of  the  discrete  PDF  vector  P#.. 
Equivalently,  we  may  write 

2  •  JSD(PRi,PR] )  =  2fT(tpfll  +  ip*2)  -  P(PHl)  -  P(Pfl2).  (3.3) 

With  the  knowledge  that  all  3D  objects  will  manifest  a  trend  similar  to  that  illustrated  in 
Figure  3.1,  we  proceed  to  determine  the  smallest  (characteristic)  resolution  5ft  of  any  given 
object  O.  This  resolution  5ft  will  also  have  a  distribution  of  its  geodesic  function  sufficiently 
close  to  that  of  the  finest  resolution.  For  a  systematic  comparison  of  the  PDFs  at  various 
resolutions,  we  first  define 

X(Ri)  =  (JSD(PR<),PRi),JSD(PRo,PRi+1 )),  i  =  1,2, . . . ,  K  —  1.  (3.4) 

We  subsequently  proceed  to  define  the  characteristic  resolution  5ft  as 

5ft  =  argmax ^variance  (X(P^))  ^ ,  i  =  1,  2, . . . ,  K  —  1.  (3.5) 

Note  that  this  resolution  evaluation  is  carried  out  during  the  learning  step  in  order  to 
alleviate  the  online  computational  burden.  To  reduce  the  computational  load  during  testing, 
we  try  to  avoid,  if  at  all  possible,  calling  upon  the  first  resolution  representation,  and 
use  instead  the  characteristic  resolution  as  the  reference.  To  that  end,  the  characteristic 
resolution  5ft  and  its  two  neighbors  along  with  the  corresponding  JSD  distances  represent 
the  reference  as  detailed  in  Figure  3.2. 


3.3  Proposed  algorithm 

With  the  tools  described  in  Section  3.2  in  hand,  we  propose  a  classification  strategy 
that  heavily  relies  on  a  training  procedure  during  which  a  class  parametrization  is  achieved 
in  order  to  yield  as  a  result,  an  object  class  membership.  This  result  is  based  on  an 
elimination  principle  which  takes  advantage  of  the  designed  parametrization.  To  further 
refine  discrimination  among  objects,  we  introduce  a  post  processing  zoom-in  procedure  on 
PDF’s  comparison  and  to  focus  on  more  detailed  dissimilarities  among  objects. 
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3.3.1  Class  resolutions 

As  noted  in  Section  3.2.1,  we  evaluate  a  characteristic  resolution  9fjW  as  defined 
in  Eq.  (3.5)  for  each  training  object  O^,  i  G  D  =  {1,  2, . . . ,  M},  from  a  class  C.  For  various 
objects  of  a  class,  we  obtain  a  class  characteristic  resolution  as  being  the  maximum  of  all 
the  characteristic  resolutions  within  that  class, 

=  max  {9?^}.  (3.6) 

i=l 

To  ensure  generalizability  to  all  training  samples,  we  define  a  second  discrimination  pa¬ 
rameter  as  being  the  transition  resolution  or  the  minimum  right  neighbor,  i.e.,  the  right 
neighbor  of  the  minimal  characteristic  resolution,  referred  to  as  Recall  that  the  mea¬ 
sure  of  variance  of  JSD  at  and  provides  the  third  parameter;  thus,  we  end  up  with 
three  parameters. 

3.3.2  Thresholding 

Upon  establishing  any  comparison  within  a  class  is  carried  out  at  by 
computing  a  pairwise  JSD  among  all  training  objects  in  C.  The  highest  distance  obtained 
from  the  latest  JSDs  is  the  sup  distance  that  no  object  within  a  class  C  would  never 
exceed. The  threshold  r  is  hence  defined  as 

r  =  max{JSD(P®,pW)},  (3.7) 

where  pW,  i  G  O,  is  the  lenient  notation  for  P^c]  the  PDF  of  the  GGF  associated  to  the 
object  Oi  from  the  class  C  and  computed  at  the  class  characteristic  resolution 
When  the  order  of  is  important  (>  1000),  the  computed  PDF  of  the  GGF  is  dependent 
of  the  sampling  rate,  i.e.,  the  number  of  bins  used  to  compute  the  histogram  of  the  GGF. 
Then,  the  threshold  r,  which  was  a  scalar  in  Eq.  (3.7),  becomes  a  vector  whose  elements 
are  thresholds  computed  at  different  sampling  rates. 

3.3.3  Zoom-in  operation 

To  detect  dissimilarities  between  two  objects  Oj  and  O&,  with  j  ^  fc,  the  corre¬ 
sponding  PDFs  are  compared  over  smaller  regions  of  their  support,  which  we  focus  on  by 
a  zoom-in  operation.  To  that  end,  denote  by  L  the  maximal  support  of  the  PDFs  of  the 
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GGFs  for  the  two  objects.  We  proceed  by  looking  at  fixed  intervals  L^,  with  Li  <  L,  over 
which  we  determine  normalized  PDFs  P^\Li)  and  P^k\Li)  for  Oj  and  Ok  respectively. 
The  JSD  between  the  two  objects  is  now  a  function  of  the  interval  Li  and  is  defined  as 

/  {Li)  =  JSD(pW(Li),  P^{Li)).  (3.8) 

Establishing  a  critical  region  over  which  two  PDFs  are  most  similar,  namely, 

L0  =  argmin  (/(£;)) .  (3.9) 

Our  interval  of  interest  is  redefined  as 

X  =  L  —  L0.  (3.10) 

The  interval  A  contains  hence  the  largest  diversity  between  Oj  and  Ok-  A  reasonable 
dissimilarity  measure  V\  2  is  the  ratio  between  the  area  on  the  object,  either  Oj  or  0&, 
corresponding  to  A  and  the  total  area  of  the  same  object  or  respectively.  Recall 

that  A  G  {[a,6]  |  0  <  a  <  6  <  1},  and  define,  for  Oj,  for  instance, 


s^{\)  =  {J2^i  U(w)eA}. 

1 

Note  that 

(3.11) 

sW=sW(  ]0,1]); 

(3.12) 

therefore,  V\  is  defined  as 

v,„o\  l(sU)W 

vx{Oj,ok)- s(j)  +  s(k)  y 

(3.13) 

In  a  more  particular  formulation  of  the  problem,  the  zoom-in  operation  introduces  a  new 

distance  measure  V ^  between  an  object  Oj  and  a  class  C  by  taking  the  second  object  Ok 

as  the  best  representative  of  the  class  C  in  the  sense  of  the  nearest  neighbor  relatively  to 

Oj  and  where  JSD  is  the  considered  distance.  Hence,  we  set  Ok  —  Oc ,  where  Oc  is  the 

object  from  class  C  of  index  l  E  Q  that  corresponds  to  min/G^{  JSD{P^\  P®)}.  It  follows 

j$n 

that 

vC{Oj)=Vx(Oj,Oc).  (3.14) 

The  two  measures  V\  and  V ^  are  dependent  of  the  fixed  length  of  the  interval  L^,  i.e, 
dependent  of  the  length  of  A.  Indeed,  the  smaller  the  length  of  A  is,  the  more  emphasis  is 
put  into  details. 

2 The  subscript  A  indicates  that  the  measure  is  achieved  for  the  fixed  interval  length  of  A. 
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3.3.4  Algorithm 

The  distinct  steps  of  the  algorithm  are  sketched  in  Figure  3.4.,  and  summarized 

next. 

1.  Define  N  object  classes  {Ci,  C2, . . . ,  Cn}. 

2.  In  the  learning  phase,  associate  to  each  class  Ci  a  corresponding  the  class 

characteristic  variance  Gq.  and  the  threshold  rci . 

3.  Sort  all  classes  in  an  increasing  resolution  order 

4.  Construct  super-classes  by  merging  classes  sharing  the  same  parameters  and 

5.  Start  from  the  lowest  resolution  and  set  l  =  1. 

6.  Compute  the  GGF  of  O  at  5?^  and  get  its  resolution  parameters. 

7.  Compare  the  resolution  parameters  of  O  with  those  of  C/.  If  similarity  is  established, 
i.e.,  decision  is  1; 

•  Termination  of  search  if  Ci  is  a  class. 


Figure  3.3:  Detection  of  the  area  of  dissimilarity  between  two  airplanes  using  the  GGF: 
Zoom-in  operation. 
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Figure  3.4:  Algorithm  of  the  discrete  classification  decision  for  3D  objects.  Resolutions 
are  simply  denoted  lRCi  in  order  not  to  overload  the  figure. 
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•  Apply  thresholding  if  Ci  is  a  super-class. 
Otherwise  the  decision  is  0  and  1  =  1  +  1. 

8.  Go  to  6. 

9.  Repeat  operations  until  a  decision  1  is  reached. 

3.4  Experimental  results 

3.4.1  Independence  of  sources  of  data 


Figure  3.5:  Illustration  of  the  characteristic  resolution  extraction  (Best  visualized  in  color). 


An  important  characteristic  of  classification  algorithm  is  its  consistent  performance 
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Table  3.1:  Characterization  parameters  for  the  class  of  felines. 


Parametrization  from  learning 

Parametrization  from  testing 

Absolute  error 

1942 

1897 

2.32  % 

% 

1501 

1420 

5.40% 

a2 

418000 

386000 

7.66% 

independently  of  what  data-base  or  measurement  system  the  data  originates  from.  So  as 
long  as  the  format  is  in  readable  form,  we  have  shown  that  our  proposed  approach  enjoys 
such  a  property.  We  made  extensive  use  of  Princeton’s  benchmark  [16]  for  all  our  learning 
and  testing.  We  have  subsequently  tested  data  from  INRIA’s  benchmark  to  validate  our 
claim.  Decision  tolerance  is  set  to  10%  as  a  maximum  absolute  error  on  each  parameter. 
The  corresponding  overall  recognition  performance,  i.e.,  number  of  correct  decisions  overall 
the  testing  experiments,  is  93.88%.  To  verify  and  illustrate  the  importance  of  the  invariance 
properties  provided  by  the  GGF,  we  applied  isometric  transforms  (rigid  and  non  rigid)  to 
the  same  dataset.  The  results  of  this  new  experiment  show  a  classification  performance 
of  97%.  The  parameters  obtained  for  the  class  of  felines,  where  a  good  representative  is 
shown  in  Figure  3.1.,  are  gathered  in  the  following  Table  3.1:  We  show,  in  Figure  3.5, 
additional  examples  of  detecting  the  characteristic  resolution  for  different  classes  of  3D 
models.  We  provide  the  comparison  of  PDFs  with  the  JSD  in  the  rightmost  column.  We, 
again,  confirm  this  trend  in  the  progression  of  the  distribution  of  the  GGF ;  indeed,  for  the 
same  color  mapping  of  the  GGF  on  the  surface  of  shapes,  we  clearly  see  a  sudden  change  in 
the  coloration  of  the  model.  This  abrupt  transition,  corresponds  each  time  to  going  below 
the  characteristic  resolution  of  a  model/class. 

3.4.2  Threshold  implementation 

Recall  the  thresholding  procedure  is  introduced  as  a  complementary  and  refined 
discrimination  measure  between  two  or  more  classes.  For  two  closely  related  classes  (e.g. 
pigs  and  dogs)  for  example,  a  resolution-based  comparison  fails  while  the  threshold-enhanced 
immediately  eliminates  the  confusion.  In  Figure  3.6.,  the  “pig”  object  is  eliminated  from  the 
class  of  dogs  because  the  maximum  pairwise  comparison  of  its  GGF  distribution  exceeds, 
over  a  region,  the  threshold  previously  computed  for  the  class  of  dogs. 
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Threshold  for  the  class  of  dogs 


Figure  3.6:  Illustration  of  the  classification  level:  Thresholding. 


3.4.3  Zoom-in  example 

The  GGF  is  sufficiently  powerful  to  detect  and  measure  dissimilarities  even  for 
relatively  close  3D  objects.  In  Figure  3.3.,  at  A  =  1/10,  we  look  at  the  area  of  highest 
difference  between  two  airplanes.  In  (a),  a  commercial  airplane  is  represented.  In  (b),  a 
more  recent  model  of  commercial  airplane  is  shown.  The  zoom-in  technique  highlights  in  (c) 
(in  white)  the  A-surface  of  difference.  This  difference  restricted  to  the  fuselage,  is  intuitively 
pleasing  as  it  confirms  our  visual  interpretation  of  a  fatter  first  airplane. 

3.5  Conclusion 

In  this  Chapter,  we  presented  a  new  3D  object  classification  strategy  based  on  four 
parameters.  We  have  shown  that  it  generally  provides  a  quick  discrete  decision  on  object 
comparison;  indeed,  once  the  learning  stage  completed,  the  classification  becomes  a  simple 
comparison  of  numbers  with  a  given  tolerance.  The  technique  is  based  on  the  distribution 
of  the  GGF  which,  in  turn,  provides  very  interesting  properties  such  as  its  invariance  to 
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all  isometric  transformations;  moreover,  a  new  and  precise  dissimilarity  measure  has  been 
introduced.  More  work  is  currently  directed  on  the  exploration  of  the  theoretical  bases  of 
the  characteristic  resolution  and  its  relation  to  topological  changes  of  surfaces. 
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Chapter  4 

3D  Shape  Partitioning 


To  simplify  the  matching  and  recognition  of  3D  objects,  we  propose  to  decompose 
a  complex  3D  shape  into  simpler  primitive  parts.  Our  proposed  partitioning  of  objects  relies 
on  their  topological  Reeb  graphs.  Taking  advantage  of  the  properties  of  Morse  theory,  we 
detect  the  critical  points  of  the  GGF.  These  points  define  the  levels  at  which  the  segmen¬ 
tation  happens.  To  preserve  the  geometry  of  objects,  we  choose  to  use  level  curves  instead 
of  intervals.  To  proceed  with  object  matching,  we  propose  a  kernel-based  technique  to 
register  Reeb  graphs.  This  optimal  positioning  of  two  Reeb  graphs  prepares  for  a  pairwise 
comparison  of  the  geometry  of  their  primitives. 

4.1  Introduction 

It  is  widely  believed  that  perception  utilizes  crucial  topological  characteristics  of 
objects.  Recent  neuro-imaging  studies,  together  with  behavioral  studies,  provide  strong 
evidence  supporting  the  notion  that  topological  properties  (such  as  a  genus,  or  connected¬ 
ness)  are  primitives  of  visual  representation  in  humans.  To  address  the  potentially  complex 
topology  of  a  3D  object  (e.g.,  holes,  breaks  in  objects),  we  call  upon  Morse  theory  [19] 
to  develop  a  robust  and  systematic  technique  which  unfolds  the  topology  of  an  object,  by 
studying  the  critical  points,  i.e.,  maximum,  minimum  or  saddle  points,  of  a  Morse  function 
defined  on  its  surface  [22];  the  GGF  in  our  case  (Section  2.3).  These  critical  points  along 
with  the  regular  points,  fully  capture  the  topology  of  a  3D  object  by  a  graph,  referred  to  as 
a  Reeb  graph.  Graph  theory  has  enjoyed  a  wide  popularity  in  problems  of  classification  and 
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recognition  [26].  Reeb  graphs  have  thus  proven  to  be  instrumental  in  globally  describing 
three-dimensional  shapes  and  in  representing  their  topology  [37].  Specifically,  the  nodes 
on  such  graphs  represent  critical  points  of  a  Morse  function  defined  on  a  3D  shape/object 
surface,  while  the  edges  capture  the  topologically  homogeneous  parts,  which  correspond  to 
“primitive”  shapes. 

Comparing  topological  graphs/skeletons  for  the  purpose  of  matching  objects  has 
been  extensively  studied.  Our  interest  is,  however,  to  explore  these  graphs  for  more  refined 
comparisons.  We  thus  propose  to  use  the  Reeb  graph  of  an  object  to  decompose  it  into 
simple  shapes  or  primitives.  Comparing  the  geometry  of  these  shapes  will  therefore  become 
an  easier  task  as  long  as  we  ensure  a  good  correspondence  between  the  primitives  of  two 
objects  to  compare. 

The  remainder  of  this  Chapter  is  organized  as  follows:  We  describe,  in  Section  4.2, 
our  original  technique  for  the  extraction  of  Reeb  graphs  using  level  curves  of  the  GGF. 
In  Section  4.2.2,  we  present  our  partitioning  technique,  and  propose,  in  Section  4.3,  a 
framework  for  object  comparison  based  on  matching  their  graphical  representation  and 
the  comparison  of  their  primitives.  Finally,  in  Section  4.4,  we  illustrate  the  results  of  our 
proposed  technique. 

4.2  Extraction  of  Reeb  graphs 


Figure  4.1:  Illustration  of  Reeb  graph  extraction:  (a)  Initial  object,  (b)  GGF  function  on 
the  surface,  (c)  Iso-geodesic  curves,  (d)  Extracted  topological  Reeb  graph.  (Best  visualized 
in  color) 
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With  the  GGF  g(-)  in  hand  (Section  2.3),  we  construct  a  topological  Reeb  graph 
for  an  object  surface  S  [1,  2].  Mathematically,  a  Reeb  graph  is  a  quotient  space  <S/~,  where 
the  equivalence  relation  is  given  by  p  ^  q  if  and  only  if  g( p)  =  g{ q)  with  p  and  q  being 
two  points  on  S  and  belonging  to  the  same  connected  component  of  g~1(g(p))-  Practically, 
we  sample  the  surface  S  of  a  3D  object  by  way  of  level  sets  of  a  Morse  function,  and  in  our 
case  thereof  a  GGF  g(-).  As  illustrated  in  Figure  4.1  (c),  the  level  sets  of  g(-)  are  closed 
curves,  which  we  refer  to  as  iso-geodesic  curves.  We  present  the  resulting  Reeb  graph  in 
Figure  4.1  (d),  where  each  closed  iso-geodesic  curve  is  replaced  by  a  node  colored  in  black. 


4.2.1  Level  set  characterization  of  a  surface 

Generically,  we  can  reconstruct  a  2-dimensional  surface  S  from  the  set  of  1- 
dimensional  iso-curves  of  the  GGF.  A  surface  S  is  then  a  disjoint  union  of  all  iso- geodesic 
sets  A{t)  for  t  G  [go  A]  C  M+,  where  A(t)  —  {v  E  S  \  g(v)  —  t},  and  go  is  the  minimal  value 

of  </(•)• 

t= 1 

5  =  |^J  A(t),  and  A(h)  D  A(t2)  =  0  if  £17^2-  (4.1) 

t=go 

For  smooth  and  compact  objects,  an  iso-geodesic  set  is  the  union  of  closed  curves  (z.e.,  iso¬ 
geodesic  curves).  The  number  of  these  distinct  curves  at  the  same  geodesic  level  is  the 
cardinality  Car{t)  of  the  corresponding  iso-geodesic  set  A(t).  In  practice  and  on  triangu¬ 
lated  meshes,  the  extraction  of  an  iso-geodesic  set  for  a  given  level,  say  A,  entails  first  finding 
all  the  faces  that  cover  this  set  of  curves.  Each  iso-geodesic  curve  is  included  in  a  connected 
set  of  covering  faces  so  that  the  resulting  curve  is  ensured  to  be  closed.  A  covering  face 
whose  vertices  are  vi,  V2  and  V3  with  g(vi)  <  g(v 2)  <  g(v3)  falls  under  one  of  the  following 
properties: 

1.  A  G  b(vi),5(v2)] 

or 


2.  A  €  b(v2),Sf(v3)] 


Two  points  p  and  q  belonging  to  the  approximated  iso-geodesic  set  M(A)  are  defined  on 
two  edges  of  the  covering  face  (Figure  4.2).  If  Property  1  is  verified,  then: 


vip 


A  -  g(y  1)  — > 

g(v2)  -  g(y  1)  VlV2’ 


26 


Figure  4.2:  Covering  face  and  iso-geodesic  curve  interpolation. 


and 


viq  = 


If  instead  Property  2  is  verified,  then: 


and 


v3p 


v3q  = 


A  -  ff(vi) 
g(v 3)  -g(v  1) 

g(y 3)  -  A 

g(v3)  -  g(v2) 

g(V3)  -  A 


V1V3. 


V3V2, 


V3V1. 


g(V3)  -  g(v l) 

Once  we  obtain  all  the  possible  couples  (p,q),  we  smoothly  interpolate,  with  a  B-spline 
for  instance,  the  so  obtained  sample  of  points  to  get  the  iso-geodesic  curve  at  the  level  A 
of  the  GGF;  moreover,  the  subsurface  supported  by  one  edge  has  exactly  one  iso-geodesic 
curve  for  every  value  between  a  and  6;  the  extremal  values  of  the  GGF  along  the  edge.  We 
call  such  subsurface  M  a  mono- cardinality  subsurface.  We  may  view  Ai  as  follows: 

M  —  | J  C'(f),  with  go  <  a  <  b  <  1, 

te[a,b] 


and  C(ti)  H  Cfa)  =  0  if  t\  ^  £2, 


(4.2) 


where  C{t)  is  just  one  iso-geodesic  curve  at  the  level  t  of  the  GGF.  In  the  example  of 
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Figure  4.3:  Illustration  of  topological  change  of  iso-geodesic  curves.  The  point  p  is  the 
bifurcation  point  at  which  this  change  occurs. 


Figure  4.3,  the  illustrated  iso-geodesic  sets  correspond  (from  top  to  bottom)  to  the  values 
0.8,  0.85  and  0.9  of  the  GGF,  respectively.  The  corresponding  cardinalities  are  Car(0.80)  = 
Car  (0.85)  =  1,  and  Car  (0.90)  =  2,  i.e.,  the  first  two  sets  consist  of  one  closed  curve,  whereas 
the  third  set  includes  two  distinct  closed  curves. 

Tracking  the  cardinalities  of  iso-geodesic  sets  helps  us  in  constructing  an  associated 
Reeb  graph.  The  continuous  and  smooth  evolution  of  the  iso-geodesic  curves  on  a  surface 
captures  its  topological  and  geometrical  description.  Any  change  in  the  topology  of  the  iso¬ 
geodesic  sets,  i.e.,  any  smooth  change  in  their  cardinalities,  hence,  determines  a  bifurcation 
point  as  shown  in  Figure  4.3.  This  transition  is  translated  on  the  object’s  topological  graph 
by  a  node  of  order  p  >  2  that  introduces  p  edges.  We  preserve  all  such  nodes  and  those 
corresponding  to  critical  points  of  the  GGF  along  with  their  edges  to  only  capture  the 
important  topological  information.  The  resulting  special  skeleton,  that  we  refer  to  as  a 
reduced  3D  Reeb  graph ,  offers  a  nice  structure  on  which  the  geometry  of  an  object  may 
easily  be  superposed.  Thus,  the  geometric  modeling  of  a  complex  shape  is  now  reduced  to 
separately  modeling  a  partial  shape  along  each  edge  of  the  reduced  Reeb  graph.  We  next 
discuss  the  details  of  modeling  the  geometric  shape  along  an  edge. 
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During  the  graph  extraction  operation,  we  replace  each  iso-geodesic  curve  Cj  by 
one  node  Nj  which  is  the  arithmetic  mean  of  the  curve;  hence,  we  end  up  with  a  point  cloud 
in  3  dimensions.  In  order  to  systematically  get  a  final  topological  representation  of  the  3D 
object,  we  define  a  canonical  relationship  ★  between  every  two  nodes  Nj  and  Nj  such  that: 

{1  if  Nj  and  Nj  are  connected, 

0  ,  otherwise. 

To  that  end,  we  give  the  following  postulates: 

•  Two  nodes  representing  the  same  iso-geodesic  level  are  disconnected.  We  write: 

i  =  j  =>  Ni*Nj  =  0.  (4.2) 

•  Two  nodes  separated  by  nodes  at  intermediate  levels  are  disconnected,  i.e., 

j  ^  i  d=  l  =>  Nj  *  Nj  =  0,  (4-2) 


where  l  is  the  sampling  step. 

Based  on  these  postulates,  we  constrain  our  search  on  nodes  representing  two  different  but 
consecutive  iso-geodesic  levels. 

Theorem  1  (Path  connectivity)  Nodes  representing  two  consecutive  iso-geodesic  curves 
are  linked  (connected)  if  and  only  if  there  exists  a  continuous  path  joining  the  two  curves 
and  lying  on  the  section  limited  by  the  same  iso-geodesic  curves. 

4.2.2  Graph  connectivity 

To  partition  a  3D  object  we  need  to  first  accurately  connect  all  the  nodes  Nj.  To 
that  end,  we  start  by  defining  orthogonal  curves. 

Definition  2  (Orthogonal  curves)  An  orthogonal  curve  on  a  surface  S  passing  through 
a  point  p  is  the  curve  of  minimal  length  linking  the  iso-geodesic  curve  containing  p  to 
another  iso-geodesic  curve. 

We  note  that  the  vector  field  generated  by  the  gradient  of  the  GGF  gives  integral  curves 
orthogonal  to  the  iso-geodesic  curves.  So,  by  definition,  an  orthogonal  curve  takes  a  point  p 
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on  a  surface  S  and  geodesically  projects  it  on  another  iso-geodesic  curve.  By  considering  an 
infinitesimal  patch  around  a  point  p  from  5,  we  approximate  the  patch  by  a  disk  on  which 
the  iso-geodesic  curve  C  becomes  a  segment  passing  through  p.  This  segment  represents 
the  direction  of  zero  variation  of  the  GGF.  We  find  that  the  projection  of  p  on  the  next 
iso-geodesic  segment  Cf  follows  the  perpendicular  to  C  on  p.  Under  the  assumption  that  all 
points  are  uniformly  distributed  on  the  surface  and  since  the  iso-geodesic  curve  represents 
the  direction  of  zero  variation  of  the  GGF,  we  conclude,  by  duality,  that  the  orthogonal 
projection  of  p  is  equivalent  to  finding  the  direction  e  of  highest  variation  of  the  GGF ;  hence, 
we  construct  an  orthogonal  curve  by  progressively  tracking,  at  a  point  level,  the  direction 
of  the  highest  variation  of  the  GGF.  We  determine  the  direction  e  that  maximizes  the 
directional  derivative  of  g  at  a  point  p;  thus,  we  extract  an  orthogonal  curve  by  progressively 
estimating  a  new  direction  e  starting  at  each  new  iso-geodesic  level1.  We  may  write, 

e  =  argmax  [Dg  •  e(p)) 

g(p  +  te)-g(p) 

—  argmax(lim - - - ).  (4.3) 

Starting  from  a  point  with  the  lowest  value  of  the  GGF,  and  by  progressively  finding  e, 
we  construct  an  orthogonal  curve  with  respect  to  iso-geodesic  curves.  By  construction, 
iso-geodesic  curves  are  transversal  to  an  orthogonal  curve.  The  orthogonal  curve  is  the 
continuous  path  that  we  were  after  in  Section  4.2.  Its  existence  between  two  consecutive 
GGF  levels  implies  the  existence  of  an  edge  connecting  the  nodes  representing  these  levels. 

The  subsurface  supported  by  one  edge  has  exactly  one  iso-geodesic  curve  for  every 
value  between  a  and  6;  the  extremal  values  of  the  GGF  along  the  edge.  We  call  such 
subsurface  M  a  mono- cardinality  subsurface.  We  may  view  M  as  follows: 

M  =  |^J  C'(t),  with  go  <  a  <  b  <  1, 

te[a,b] 

and  C(t\)  fl  Cfo)  =  0  if  t\^t^  (4.4) 

where  C(t)  is  just  one  iso-geodesic  curve  at  the  level  t  of  the  GGF.  In  the  example  of 
Figure  4.3,  the  illustrated  iso-geodesic  sets  correspond  (from  top  to  bottom)  to  the  values 
0.8,  0.85  and  0.9  of  the  GGF,  respectively.  The  corresponding  cardinalities  are  Car(0.80)  = 


2We  note  that  we  only  need  a  starting  point  to  extract  an  orthogonal  curve. 
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Starting  point  p  at  the 


Figure  4.4:  Orthogonal  curve  extraction  using  infinitesimal  patches  around  each  point  at 
increasing  values  of  the  GGF. 


Car  (0.85)  =  1,  and  Car  (0.90)  =  2,  i.e.,  the  first  two  sets  consist  of  one  closed  curve,  whereas 
the  third  set  includes  two  distinct  closed  curves. 

Tracking  the  cardinalities  of  iso-geodesic  sets  helps  us  in  constructing  an  associated 
Reeb  graph.  The  continuous  and  smooth  evolution  of  the  iso-geodesic  curves  on  a  surface 
captures  its  topological  and  geometrical  description.  Any  change  in  the  topology  of  the  iso¬ 
geodesic  sets,  i.e.,  any  smooth  change  in  their  cardinalities,  hence,  determines  a  bifurcation 
point  as  shown  in  Figure  4.3.  This  transition  is  translated  on  the  object’s  topological  graph 
by  a  node  of  order  77  >  2  that  introduces  77  edges.  We  preserve  all  such  nodes  and  those 
corresponding  to  critical  points  of  the  GGF  along  with  their  edges  to  only  capture  the 
important  topological  information.  The  resulting  special  skeleton,  that  we  refer  to  as  a 
reduced  3D  Reeb  graph ,  is  the  main  support  for  a  meaningful  segmentation  of  a  3D  object. 
Indeed,  the  only  remaining  task  is  to  cluster  all  the  points  on  the  surface  S  that  are  delimited 
by  two  iso-geodesic  curves  represented  by  two  connected  nodes  from  the  reduced  graph. 
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4.2.3  Efficient  sampling 

A  very  redundant  and  important  problem  that  is  stressed  when  building  Reeb 
graphs  is  how  to  solve  this  blind  setting  and  define  the  effective  sampling  rate  K  that 
ensures  the  extraction  of,  and  only  of,  the  information  needed  to  represent  and  reconstruct 
a  given  shape  [1,  28].  The  sampling  rate  K  is  usually  empirically  chosen,  and  the  resulting 
sampling  step  is: 

_  maxves  00))  -  minv&s  (</(v)) 

K  '  1 

The  use  of  curvature  as  a  criterion  is  certainly  the  first  thing  that  comes  to  one’s  mind. 
Indeed,  if  we  think  of  planar  shapes,  we  know  that  more  points  (vertices)  are  needed  to 
represent  high  curvatures  in  opposition  to  straight  line  segments  whose  curvature  is  zero 
and  for  which  the  two  extremal  points  are  sufficient  to  represent  and  reconstruct  the  full 
shape.  Our  intuition  may  further  suggest  to  take  the  highest  curvature  on  the  curve/surface 
to  be  the  criterion.  This  would  be  the  solution  if  the  highest  curvature  happened  to  be  a 
dominant  feature  on  the  whole  shape.  In  many  cases,  this  is  not  true.  Instead  of  basing  our 
choice  on  a  local  property,  we  therefore  rely  on  the  global  perception  of  shapes’  curvatures. 

The  characteristic  resolution  5ft,  defined  in  Eq.  (3.3),  is  a  feature  that  is  directly 
related  to  the  curvature  of  a  shape.  As  we  reduce  5ft,  we  act  on  the  original  shape  by 
smoothing  it.  This  machinery  is  nothing  but  choosing  a  tolerance  r.  “Tolerance”  is  the 
common  name  given  for  the  maximal  edge  length  between  adjacent  points  on  a  surface 
S.  In  fact,  the  raw  data  is  reduced  to  its  characteristic  resolution  5ft  and  every  vertex 
becomes  now  of  importance  to  obtain  an  accurate  representation  for  S.  We  recall  that  any 
topological  transformation  that  happens  on  S  and  that  we  are  to  detect  is  translated  on  its 
one  dimensional  Morse  function  g(-);  hence,  we  base  the  sampling  rule  of  a  3D  object  on  the 
sampling  of  g(-)  or  of  a  functional  of  g(-).  We  view  g(-)  as  a  continuous  random  variable, 
or  practically  as  a  discrete  random  variable  X  (g  =  X),  s.t.,  X  :  V  — >  M  is  defined  on  the 
probability  space  (V,A,  Psr),  where  V  is  the  set  of  5ft  vertices  on  the  triangle  mesh  (V,  .F) 
approximating  the  surface  S  (Section  2.1).  A  is  a  a  —  Algebra  composed  of  all  the  subsets 
of  V.  Note  that  we  interchangeably  use  V  as  the  set  or  the  matrix  of  vertices. 

Similarly  to  Eq.  (4.5),  we  apply  the  sampling  along  the  axis  of  variation  of  X. 
In  order  to  achieve  an  efficient  sampling  step  /,  we  need  to  define  the  largest  step  l  that 
will  cause  no  breakdowns  on  the  mesh  (V,F).  Let  71  =  g~x  (g( p))  and  72  =  g~l  (ff(q)) 
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be  two  distinct  and  consecutive  connected  components  (i.e.,  two  consecutive  iso-geodesic 
curves).  With  a  resolution  3?  and  a  tolerance  r  associated  to  it,  the  longest  acceptable 
orthogonal  geodesic  distance  between  71  and  72  is  equal  to  r.  The  remaining  question  is  to 
find  l  —  \g(p )  —  g( q)|  =  |Ai  —  A2 1  -  If  p  and  q  are  connected  then  d  (p,  q)  =  r,  and  we  have 
for  v  G  V: 

d(p,  v)  >  d(p,  v)  =4>  d(p,  v)  =  d(q,  v)  +  r,  (4.6) 


Similarly, 

d(p,  v)  <  d(p,  v)  =4>  d(q,  v)  =  d(p,  v)  +  r; 
Therefore,  if  there  are  n  elements  from  V  verifying  (4.6),  then  we  find: 

(R-2(n-l)) 


Ai  -  A2  = 


maxvGV  Ewev  d  (v? 


•  T. 


(4.7) 


(4.8) 


We  may  also  define  an  upper  bound  for  /,  independent  of  the  variable  n  and  equal  to 

(K-2) 

maXvGV  Ewev  d(v>w) 


4.3  Reeb  graph  matching 


Reeb  graphs  being  topological  descriptors,  we  require  the  connectivity  information 
to  define  a  distance  between  them.  We  actually  use  reduced  Reeb  graphs  where  the  node 
order  is  never  2.  We  herein  define  a  distance  measure  that  is  in  the  same  spirit  as  the 
edit-distance  for  shape  matching  as  defined  in  [4].  We  compute  in  (4.9)  the  similarity  score 
§(•,  •)  between  two  reduced  Reeb  graphs  T i  and  T2  with  N±  and  7V2  nodes,  respectively,  and 
N\  ^  7V2.  We  denote  by  {v^}^  the  nodes  of  Ti,  and  {w i}^1  the  nodes  of  T2.  We  define 
the  similarity  between  these  two  graphs  as  follows: 

2  Nl 

S(ri,  r2)  =  = may  {^(vi,  wi)<Jp)(wi)},  (4.9) 

2=1  J 


with 


and 


(wi) 


1  iff  | r/v,  -  rjwj  I  <  2; 
0  otherwise, 


/?(vj,  wj) 


1 

2n<j 


exp 


b(vi)  -^(wj)| 

2cj2 


(4.10) 
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where  a  is  a  strictly  positive  constant  of  our  choosing  and  rjWi  and  r]^Vj  are  the  orders  of 
nodes  and  Wj,  respectively.  The  function  S^Tk\‘)  plays  the  role  of  a  registration  operator 
based  on  the  order  of  graph  nodes.  We  also  choose  a  kernel  function  /?(-,  •)  that  is  strictly 
decreasing,  hence  the  choice  of  the  Gaussian  kernel  in  Eq.  (4.10).  The  intuition  behind 
choosing  the  factor  (3  is  the  assignment  of  a  high  similarity  weight  to  registered  nodes  if 
they  correspond  to  GGF  levels  that  are  close. 

4.4  Illustration 

In  order  to  illustrate  the  proposed  shape  partitioning  technique,  we  use  the  four 
objects  in  Figure  4.5  (a)  [15].  We  start  by  computing  the  GGF  corresponding  to  each  one 
of  these  shapes  and  extract  the  corresponding  topological  Reeb  graphs  (Figure  4.5  (b)).  In 
Figure  4.5  (c),  we  present  the  segmented  objects.  These  have  the  same  colors  for  primitives 
that  we  matched  using  Eq.(4.9). 

4.5  Conclusion 

We  proposed  a  technique  for  an  automatic  partitioning  of  3D  objects  into  prim¬ 
itives  along  with  a  distance  to  ensure  a  good  correspondence  between  them.  The  next 
Chapter  shows  how  this  segmentation  may  be  used  to  simplify  recognition  problems  and 
for  partial  matching  of  shapes. 
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Figure  4.5: 


Overview  of  the  partitioning  technique  (Best  visualized  in  color). 
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Chapter  5 

Squigraphs  for  Geometric 
Modeling 


We  propose  to  superpose  global  topological  and  local  geometric  3D  shape  descrip¬ 
tors  in  order  to  define  one  compact  and  discriminative  representation  for  a  3D  object.  While 
a  number  of  available  3D  shape  modeling  techniques  yield  satisfactory  object  classification 
rates,  there  is  still  a  need  for  a  refined  and  efficient  identification/recognition  of  objects 
among  the  same  class.  In  this  Chapter,  we  add  detailed  geometric  information  to  Reeb 
graphs,  by  tracking  the  evolution  of  Morse  function’s  level  curves  along  each  primitive  re¬ 
sulting  from  partitioning  complex  shapes.  We  then  embed  the  manifold  of  these  curves  into 
M3,  and  obtain  a  single  curve.  We  thus  build  new  graphs  rich  in  topological  and  geometric 
information  that  we  refer  to  as  squigraphs.  Our  experiments  show  that  squigraphs  are  more 
general  than  existing  techniques.  They  achieve  similar  classification  rates  to  those  achieved 
by  classical  shape  descriptors.  Their  performance,  however,  becomes  clearly  superior  when 
finer  classification  and  identification  operations  are  targeted.  Indeed,  while  other  techniques 
see  their  performances  dropping,  squigraphs  maintain  a  performance  rate  of  the  order  of 
97%. 


5.1  Introduction 

Numerous  3D  representation  approaches  have  recently  been  proposed,  as  overviewed 
in  [29]  and  [30].  In  the  present  work,  we  propose  a  3D  object  representation  model  that 


36 


offers  different  levels  of  discrimination.  This  model  invokes  an  object’s  geometric  and  topo¬ 
logical  information.  This  rather  recent  technique,  henceforth  referred  to  as  topo- geometric 
modeling  of  3D  objects,  was  first  exploited  in  [31],  where  Meng  Yu  et  al.  simultaneously 
define  topological  and  geometrical  feature  maps.  They  show  that  these  feature  maps  are 
invariant  to  all  affine  transforms.  This  invariance  includes,  by  definition,  non  uniform  scal¬ 
ing  transforms.  This,  in  turn,  implies  an  inaccurate  dissimilarity  measure  between  the 
geometry  of  shapes,  as  the  geometry  of  shapes  is  only  invariant  to  Euclidean  transforms. 
Tung  and  Schmitt  [2]  and  Baloch  et  al  [28]  present  a  different  theme  by  first  representing 
the  topology  of  a  3D  shape,  and  subsequently  enhancing  it  with  its  geometrical  represen¬ 
tation.  The  advantage  of  separating  the  modeling  into  two  distinct  steps  is  twofold:  First 
it  provides  two  levels  of  discrimination;  a)  a  coarse  level  with  a  simple  topological  skele¬ 
ton  called  a  Reeb  graph  [22,  32]  (Section  4.2),  and  b)  a  fine  level  with  geometric  weights 
assigned  to  each  edge  or  node  of  the  previously  extracted  topological  graph.  Another  ad¬ 
vantage  of  a  topological  graph  representation  of  an  object  is  its  ability  of  matching  objects 
by  parts,  thereby  enabling  a  transition  from  a  global  to  a  localized  correspondence  between 
shapes  (e.g.,  mechanical  parts,  manufactured  solids).  In  representing  an  object,  one  ideally 
avoids  choosing  a  reference  point  extrinsic  to  its  surface  [1,  33,  34,  35,  36].  To  that  end, 
we  address  this  limitation,  evident  in  [28,  31],  by  adopting  an  integrated  geodesic  distance 
function  as  defined  in  [1].  Such  a  function  ensures  the  invariance  of  the  Reeb  graph  of  an 
object  subjected  to  any  isometric  transform.  As  first  proposed  in  [1],  this  function  was  lim¬ 
ited  to  only  providing  topological  information  of  an  object  surface.  Tung  and  Schmitt  [2], 
in  an  attempt  to  mitigate  this  limitation,  proposed  to  revert  back  to  a  Euclidean  reference 
frame  to  compute  differential  geometric  features  which  they  use  as  weight  attributes  on  the 
Reeb  graph. 

In  this  Chapter,  we  propose  a  novel  alternative  technique  that  is  efficient,  theoretically 
sound,  practically  complete,  and  computationally  simple.  Our  proposed  approach  is  rooted 
in  modeling  the  geometry  along  edges  of  the  object  representative  graph.  We  subsequently 
define  a  new  topo-geometric  skeleton  that  we  call  squigraph.  The  key  idea  for  the  geometric 
modeling  of  shapes  is  to  embed  in  Euclidean  space  a  manifold  of  new  characteristic  curves 
(we  refer  to  as  iso-geodesic  curves)  of  an  object  surface.  The  resulting  embedding  space,  as 
we  show  below,  is  R3,  and  the  curve  manifold  is  1-dimensional,  i.e.,  a  space  curve  in  R3. 
As  a  result,  all  topological  and  geometric  information  we  exploit  remains  intrinsic  to  the 
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modeled  object,  hence  preserving  all  desired  invariance  properties  alluded  to  earlier. 

5.2  Representation  of  a  surface  geometry 

We  note  that,  just  as  in  [2,  28],  our  3D  object  modeling  approach  starts  with  a 
topological  analysis  followed  by  a  geometrical  analysis  step.  In  our  present  work,  however, 
we  strive  to  achieve  a  topo-geometric  representation  that  is  both  discriminative  and  effi¬ 
cient.  This  is  why  we  proceed  by  economically  exploiting,  for  the  geometric  modeling,  the 
same  entities  used  earlier  for  the  topological  modeling  with  Reeb  graphs.  The  entities  in 
question  are  the  level  sets  of  the  GGF,  i.e.,  the  iso- geodesic  sets  defined  on  the  surface  of 
an  object.  Using  iso-geodesic  sets  allows  us  to  further  extend  the  invariance  properties  of 
the  GGF  to  the  geometric  phase. 

In  addition,  unlike  the  usual  global  shape  descriptors,  a  more  complete  geometric  repre¬ 
sentation  is  one  which  provides  a  local  description  taking  into  account  the  spatial  location 
of  points  or  vertices  [31].  The  iso-geodesic  curves  used  in  the  topological  representation 
consequently  appear,  once  again,  to  be  able  to  offer  perfect  local  geometric  descriptors.  We 
therefore  take  advantage  of  these  already  extracted  entities,  and  use  them  for  the  second 
phase  of  our  modeling  as  described  in  what  follows. 

Broadly,  our  strategy  is  as  follows:  once  we  extract  iso-levels  of  the  GGF,  we  ob¬ 
tain  a  set  of  closed  curves  (namely  the  iso-geodesic  curves)  along  each  edge  of  a  topological 
graph  of  an  object.  We  use  these  characteristic  curves  to  define  a  geometric  model.  We 
subsequently  assign  the  models  to  the  corresponding  edges.  We  in  turn  exploit  these  geo¬ 
metric  models  to  further  compare  every  two  edges  that  have  been  matched  as  part  of  two 
topologically  similar  graphs.  In  Figure  4.5  (c),  we  show  all  matched  edges  in  the  same  color. 
We  compactly  represent  the  final  topo-geometric  model  through  a  new  spatial  graph  that 
we  call  squigraph  (Figure  5.2  (d)).  A  squigraph  graph  differs  from  a  classical  Reeb  graph  by 
its  “squiggling”  curves  that  replace  the  standard  straight  edges  and  encode  the  geometric 
information  about  the  shapes. 
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5.3  Whitney  theorem  for  modeling 

It  is  important  to  note  that  a  geometric  shape  along  an  edge  corresponds  to  what 
we  described,  in  Section  4.2.1,  as  a  mono-cardinality  subsurface.  So  all  the  partial  shapes 
along  Reeb  graph  edges  are  simply  equivalent  to  generalized  cylinders.  In  fact,  this  philos¬ 
ophy  of  decomposing  a  complex  shape  into  a  set  of  cylinders  goes  back  to  the  very  early 
work  on  tubular  representation  of  3D  shapes. 

Equation  (4.4)  shows  that  a  mono-cardinality  subsurface  Ai  is  formed  by  the  disjoint  union 
of  an  infinite  number  of  iso-geodesic  curves  C(t),  monotonically  and  continuously  parameter¬ 
ized  by  t  E  [a,  b]  C  R+,  where  [a,  b\  is  the  range  of  the  GGF  restricted  to  Ai  (Figure  5.1  (a)). 
It  therefore  naturally  follows  that  At  is  a  smooth  path  in  high  dimensional  space  whose 
point  elements  are  the  iso-geodesic  curves  C(t),  t  E  [a,  6];  hence,  At  constitutes  a  smooth 
1-dimensional  manifold  on  the  topological  space  of  closed  iso-geodesic  curves  on  which  the 
restriction  of  the  GGF  is  a  homeomorphism  (Figure  5.1  (c)).  In  what  follows,  we  inter¬ 
changeably  use  M  to  refer  to  the  mono-cardinality  subsurface  embedded  in  R3,  as  well  as 
to  the  corresponding  curve  embedded  in  the  higher  dimensional  space. 

Getting  back  to  our  objective  of  simply  modeling  each  mono-cardinality  subsurface 
At,  we  first  invoke  the  Whitney  embedding  theorem  as  stated  below, 

Theorem  2  (Whitney  embedding  theorem)  [38]  Let  7  be  an  n-dimensional  compact 
Hausdorff  Cr  manifold ,  2  <  r  <  00.  Then,  there  is  a  Cr  embedding  of  7  in  R2n+1. 

In  the  present  problem,  the  manifold  7  is  nothing  but  At,  for  which  the  dimension  n  is 
equal  to  one;  hence,  the  stated  theorem  asserts  the  existence  of  an  embedding  of  At  in  R3. 
In  fact,  such  embeddings  have  full  measure  in  the  set  of  maps  or  even  projections  into  R3, 
so  we  can  assume  that  a  random  projection  works.  Thus  we  can  safely  use  this  embedding 
for  our  geometric  modeling  since,  by  definition,  an  embedding  preserves  the  geometry  and 
does  not  introduce  any  new  intersections  in  the  new  space.  In  practice,  we  approximate  At 
with  a  finite  number  N  of  iso-geodesic  curves  as  illustrated  in  Figure  5.1  (b).  On  each  curve, 
we  take  M  uniformly  spaced  points  wl-  defined  by  their  Euclidean  coordinates  (x1-,  y1-,  Zj), 
i  =  1,2,...,  2V,  and  j  =  1,2,...,  M.  Each  curve  is  now  a  point  in  R3M+1  represented  by 
the  column  vector  V^,  i  =  1,  2, . . . ,  N  (Figure  5.1  (c)),  such  that, 


T 


V  = 

vf 

T 

vif 

T 

1 

with  Vj  =  [x},y),z]\T  . 

(5.1) 

. 
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Figure  5.1:  Geometric  modeling  of  a  mono-cardinality  subsurface  Ai:  (a)  The  GGF  on 
Ai  is  strictly  monotonous  and  takes  its  values  in  the  interval  [a,  b\.  As  a  result,  the  two 
bounding  iso-geodesic  curves  are  C(a)  and  C(b).  (b)  Discretized  version  of  Ai.  (c)  Path 
created  by  Ai  in  high  dimensional  space,  (d)  Final  modeling  curve  in  3D  space. 
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Figure  5.2:  Object  representation  using  a  squigraph.  (a)  Original  shape,  (b)  Extracted  Reeb 
graph,  (c)  Partitioning  into  primitives,  (d)  Squigraph  for  a  topo-geometric  representation. 


The  sample  set  of  iso-geodesic  curves  on  At  is  a  matrix  V  of  dimension  3 M  x  N  where 
V  =  [Vi  •  •  •  Vjv]-  Applying  the  Whitney  theorem  on  V  reduces  the  embedding  problem  to 
a  simple  linear  formulation; 

W  =  PTV,  (5.2) 

where  W  is  a  set  of  N  points  in  R3  resulting  from  the  projection  of  V  via  P  into  R3.  In 
other  words,  as  illustrated  in  Figure  5.1  (d),  the  manifold  At  is  reduced  to  a  space  curve  7 
whose  sample  points  are  represented  by  the  columns  of  the  matrix  W.  In  Figure  5.2,  we 
illustrate  the  idea  of  a  topo-geometric  modeling  via  the  space  modeling  curves  that  we  just 
defined;  hence,  we  intrinsically  enhance  the  typical  3D  Reeb  graph  of  Figure  5.2  (b).  To 
that  end,  we  assign  one  modeling  curve  to  each  edge  of  the  Reeb  graph.  By  so  doing,  we 
may  view  the  final  representation  as  a  new  kind  of  graph,  that  we  refer  to  as  a  squigraph , 
as  shown  in  Figure  5.2  (d).  The  Whitney  embedding  theorem  guarantees  that  almost  any 
projection  of  At  is  an  embedding.  In  [39,  40],  the  notion  of  good  Whitney  embedding  is 
introduced  to  identify  a  class  of  projections.  For  many  applications,  such  as  recognition, 
a  unique  representation  of  an  object  is  required.  As  a  result,  we  apply  the  notion  of 
optimal  embedding  with  the  performance  criterion  presented  in  [40]  and  defined  in  (5.4) 
being  maximized. 
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5.4  Optimal  projection 

An  embedding  gives  us  a  one-to-one  mapping  between  the  path  in  high  dimensional 
space  M  C  R3M,  and  the  path  7  C  R3,  so  that  no  two  points  in  A4  will  collapse  as  a  result 
of  the  mapping  in  (5.2).  An  implementation  of  this  principle  is  presented  in  [40]  and  referred 
to  as  a  secant  based  method.  As  just  explained,  no  two  points  in  R3M  are  mapped  to  the 


Figure  5.3:  Asymptotic  convergence  of  LTMADS  applied  to  3D  shape  dataset. 

same  point  in  R3.  This  means  that  the  projecting  vector  is  linearly  independent  from  any 
possible  secant1  in  R3M.  We  define  the  set  of  all  possible  unit  secant  vectors  from  the  initial 


1  Secant  or  secant  line  is  a  line  that  intersects  two  points  from  a  curve. 
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data  V  as  follows: 


IV  _  vV  (*’•*)  G  i1’---,  W  and  i±j 
I  *  i  *  j  II 


=  {¥*,  i  =  1, . . . ,  L}  , 


N\ 


with  L  =  (iV_2)!2!- 

The  performance  of  the  projection  operator  P  is  then  reflected  by: 


(5.3) 


e  =  .jnin^  [||  P T^i  ||2]  .  (5.4) 

The  closer  e  is  to  1,  the  better  spread  out  the  projected  points  in  R3  are,  and  the  better 
the  choice  of  P  is;  hence,  finding  a  good  estimate  of  the  optimal  projection  is  reduced  to 
the  following  minimization  problem: 


P 


arg  min 


min 


arg  minP(P). 


P  h 


(5.5) 


Since  the  function  F( P)  is  non  differentiable,  we  need  a  direct  search  algorithm  that  does 
not  require  any  derivative  of  the  function  to  minimize  it.  As  presented  in  [41],  the  Lower 
Triangular  Mesh  Adapted  Direct  Search  algorithm  (LTMADS)  has  been  adapted  to  Rie- 
mannian  manifolds  and  tested  on  an  example  of  a  Whitney  embedding.  The  solution  given 
by  LTMADS  is  highly  dependent  on  the  initial  data.  On  our  3D  objects,  the  performance 
of  LTMADS  converges  asymptotically  to  1  (Figure  5.3).  We  define  a  stopping  criterion  for 
the  LTMADS  algorithm  adapted  to  each  data  by  considering  the  slope  Ae  =  \ei  —  e^_ i|, 
where  C{  and  \  respectively  represent  the  performances  of  the  optimization  algorithm  at 
iterations  i  and  i  —  1. 


5.5  Space  marking 

It  is  important  to  ensure  that  all  iso-geodesic  curves’  discrete  representations  are 
uniform  and  consistent  within  a  coordinate  space.  Visually,  we  may  equate  this  task  to 
marking  the  surface  Ml.  By  marking,  we  mean  drawing  a  well  defined  reference  line  on  the 
surface  Ml.  In  Figure  5.1  (a),  this  reference  line  is  the  one  curve  going  through  the  two  red 
points  corresponding  to  levels  a  and  b.  This  special  line  enables  us  to  construct  the  vectors 
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in  an  accurate  and  consistent  manner.  To  that  end,  we  use  the  orthogonal  curve  defined 
in  Eq.  (4.3). 

We  note  that  we  only  need  a  starting  point  to  extract  an  orthogonal  curve.  This 
point,  however,  is  not  unique  as  it  is  only  required  to  be  a  point  from  a  bounding  iso¬ 
geodesic  curve  ( C(a )  or  C(b)  in  Figure  5.1).  In  fact,  thanks  to  the  construction  proposed 
in  Eq.  (5.1),  choosing  a  different  starting  point  will  only  result  in  a  mere  rotation  of  the 
reference  axes  in  R3M.  Besides,  because  in  the  present  geometric  modeling  we  require  a  full 
invariance  to  Euclidean  transforms,  this  rotation  will  have  no  effect  on  the  final  modeling 
curve  (7  in  Figure  5.1  (d)). 
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Chapter  6 

Correspondence  Preserving 
Similarity  Invariant  for  Space 
Curves 


We  present  a  new  similarity  invariant  signature  for  space  curves.  This  signature 
is  based  on  the  information  contained  in  the  turning  angles  of  both  the  tangent  and  the 
binormal  vectors  at  each  point  on  the  curve.  For  an  accurate  comparison  of  these  signatures, 
we  define  a  Riemannian  metric  on  the  space  of  the  invariant.  We  show  through  relevant 
examples  that,  unlike  classical  invariants,  the  one  we  define  in  this  paper  enjoys  multiple 
important  properties  at  the  same  time,  namely,  a  high  discrimination  level,  independence 
of  any  reference  point,  uniqueness,  as  well  as  a  good  preservation  of  the  correspondence 
between  curves;  moreover,  we  show  how  to  use  the  proposed  signatures  for  both  full  and 
partial  matchings  of  3D  objects. 

6.1  Introduction 

Multiple  3D  modeling  methods  use  spatial  curves  for  recognition  and  matching  of 
objects.  Spatial  curves  are  exploited  in  different  configurations.  They  may,  for  instance,  be 
extracted  as  contours  of  landmark  surfaces  [7],  as  level  curves  of  a  Morse  function  [8],  or 
also  as  elements  of  curved  skeletons  [9,  10].  All  these  techniques,  despite  their  differences, 
agree  in  relying  on  curves’  properties  in  solving  computer  vision  problems.  This  common 
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approach  is  motivated  by  the  fact  that  curves  in  3D  are  fairly  well  known  geometric  entities; 
moreover,  under  some  conditions,  they  can  accurately  describe  the  overall  geometry  of  an 
object  in  3D  space  [7].  Translating  the  constraints  of  3D  shape  representation  techniques  to 
curves  reduces  the  level  of  difficulty  associated  with  the  representation  problem  and  makes 
it  more  tractable. 

Pose  invariance  of  surfaces  is  a  common  requirement  in  object  modeling.  It  is  also  a  good 
illustration  of  the  simplification  from  surfaces  to  curves.  Indeed,  an  effective  and  economical 
solution  for  curves  pose  invariance  may  be  provided  through  Euclidean/similarity  invariants 
or  invariant  signatures  [42].  Besides  this  pose  invariance  property,  additional  constraints  are 
imposed  on  3D  curves  as  a  direct  result  of  the  nature  of  3D  shape  recognition  applications. 
These  constraints  may  be  summarized  as  follows:  invariance  to  a  group  of  transforms, 
uniqueness,  local  characterization  (local  support),  ability  to  determine  shape  properties 
such  as  symmetries  and  part  correspondences.  To  the  best  of  our  knowledge,  none  of  the 
available  references  seem  to  gather  all  these  properties  at  once.  The  most  complete  work 
is  the  one  of  Mokhtarian  and  Bober  [7],  as  they  succeed  in  citing  and  addressing  all  the 
necessary  properties;  however,  to  provide  an  invariance  to  scaling  transforms,  the  authors 
use  a  multi-resolutional  procedure.  In  the  present  work,  we  provide  the  advantage  of  having 
an  invariant  that  is,  by  definition,  i.e.,  without  additional  steps,  fully  invariant  to  all  sim¬ 
ilarity  transforms.  The  key  contribution  of  our  work  resides  in  using  turning  angles  based 
on  curvature  and  torsion  instead  of  using  curvature  and  torsion  directly.  We  further  ensure 
a  natural  registration  of  all  the  invariants  on  one  curved  space  which  leads  to  defining  an 
accurate  and  computationally  easy  metric  for  curves  comparison.  Indeed,  we  show  that 
while  torsion  and  curvature  are  clearly  variant  with  scaling,  turning  angles  are  not.  The 
first  inspiration  of  our  work  comes  directly  from  [43],  where  an  invariant  for  planar  curves 
is  defined.  This  invariant  has  the  particularity  to  be  an  information  theoretical  measure 
of  local  geometric  properties  of  curves.  Moreover,  this  invariant  comes  as  a  proof  for  long¬ 
time  psychological  assumptions  on  mental  shape  perception.  The  operation  of  comparing 
invariants,  although  often  overlooked,  is  crucial  in  assessing  the  properties  of  the  invariant, 
and  accurately  achieving  recognition  operations.  This  is  why,  in  the  present  work,  we  com¬ 
plete  the  proposed  invariant  signature  by  defining  its  Riemannian  space  equipped  with  an 
intrinsic  measure. 

The  present  Chapter  is  organized  as  follows:  In  Section  6.2,  we  explain  the  mo- 
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tivation  behind  this  work.  In  Section  6.3,  we  review  important  geometric  and  information 
theoretical  notions  used  throughout  this  paper.  We  briefly  cover  the  work  of  Feldman  and 
Singh  [43]  in  Section  6.4,  as  it  was  the  precursor  and  the  inspiration  for  this  effort.  In 
Section  6.5,  we  introduce  our  new  similarity  invariant  signatures  for  space  curves  and  define 
a  Riemannian  metric  for  their  accurate  comparison  in  Section  6.6.  We  prove  and  illustrate 
the  property  of  correspondence  preservation  in  Section  6.7.  Finally,  in  Section  6.8  we  show 
direct  applications  of  the  proposed  signatures  in  3D  object  matching  and  comparison. 

6.2  Motivation  and  related  work 

We  say  that  two  space  curves  have  the  same  shape  if  there  exists  a  rigid  trans¬ 
formation  (rotation,  translation)  or  a  uniform  scaling  or  a  combination  thereof  that  will 
cause  one  of  these  curves  to  completely  overlap  the  other.  In  addition,  we  talk  about  curves 
that  are  partial  matches.  This  practically  translates  into  shapes  subjected  to  occlusions, 
or  problems  of  reconstruction  from  multiple  parts.  In  what  follows  we  provide  a  list  of 
properties  that  models  of  space  curves  are  to  verify; 

1.  Invariance:  In  this  work,  we  view  a  shape  as  a  geometrical  entity  that  is  independent  of 
the  coordinate  space.  It  is  therefore  necessary  to  have  a  representation  that  eliminates 
the  effects  of  rigid  motions  and  uniform  scalings.  In  other  words,  we  request  our 
modeling  to  be  invariant  to  the  group  of  similarity  transforms. 

2.  Uniqueness:  Curves  with  the  same  representation  are  considered  to  be  curves  with 
the  same  shape.  If  curves’  models  are  different,  then  these  curves  are  geometrically 
different.  This  requirement  implies  using  a  complete  system  of  invariant  functions 
whose  combination  will  define  the  final  model. 

3.  Local  support:  This  property  is  related  to  the  partial  matching  problem.  In  this  work 
our  objective  is  to  respect  the  correspondences  between  parts  of  curves.  This  is  why 
we  highlight  the  importance  of  having  a  locally  supported  representation. 

4.  Independence  of  parametrization  and  of  the  starting  point:  This  property  will  show 
to  be  important  in  defining  the  appropriate  metric  to  compare  the  final  models  for 


space  curves. 
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Figure  6.1:  Frenet  Frame. 


5.  Determination  of  shape  properties:  It  is  sometimes  of  interest  to  be  able  to  determine 
some  shape  properties  such  as  symmetries  or  periodicities. 

As  we  progress  in  this  paper,  we  will  relate  our  results  to  the  motivations  listed  above.  We 
start  this  work  by  identifying  appropriate  geometric  features  to  use  in  defining  our  model  for 
space  curves.  We  propose  to  rely  on  the  curvature  and  the  torsion  functions  of  a  space  curve 
since  these  two  uniquely  specify  a  space  curve  up  to  a  rotation  and  a  translation;  however, 
they  are  not  invariant  to  scaling  which  leads  us  to  using  the  corresponding  turning  angles, 
instead. 


6.3  Background  and  formulation 

6.3.1  Turning  angles 

A  space  curve  is  uniquely  determined,  up  to  a  Euclidean  transform,  by  its  curvature 
function  and  torsion  function  r(t),  both  continuous  functions  of  the  parameter  £; 
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hence,  we  naturally  use  these  measurements  to  define  an  adequate  invariant  signature  curve; 
however,  since  we  target  the  group  of  similarity  transforms,  and  knowing  that  curvature  and 
torsion  are  not  scale  invariant,  we  use  turning  angles  as  the  geometric  features  describing 
space  curves  [43].  In  what  follows,  we  show  how  it  is  possible  to  relate  curvature  and  torsion 
as  linear  functions  of  turning  angles.  Using  the  Frenet-Serret  formulae  (Figure  6.1), 


dT 

dt 

=  kN, 

(6.1) 

dN 

dt 

=  — kT  +  rB, 

(6.2) 

dB 

dt 

=  -rN, 

(6.3) 

where  T,  N  and  B  are  the  tangent,  the  normal  and  the  binormal  vectors,  respectively,  we 
define  two  turning  angles  olt  and  <a#.  olt  is  the  change  in  the  direction  of  T,  and  olb  is 
the  change  in  the  direction  of  B  such  that:  ax(t)  «  n(t)  •  dt  and  as(t)  ~  r(t)  •  dt.  By 
choosing  to  use  curvatures  as  main  descriptors,  we  systematically  verify  properties  1)  and 
2)  in  the  list  of  motivations,  as  well  as  the  independence  of  parametrization.  We  also  start 
touching  upon  properties  3)  and  4). 

6.3.2  Shannon  surprisal 

We  extract  the  information  contained  in  the  proposed  turning  angles,  i.e.,  our 
chosen  geometric  measures,  by  using  a  notion  in  information  theory  known  as  “Shannon 
surprisal”  [44,  43].  We  assume  that  a(t)  follows  a  von  Mises  distribution  with  a  zero  mean 
and  a  spread  parameter  b  equal  to  1.  We  then  define  the  probability  density  function  of 
a(t)  as  follows: 

fa  (a(t))  =  A  exp  (cos (a(£))) ,  V  £,  (6.4) 

with  a  E  [—7 r,  +7r]  and  A  —  ^BesseMb  b)  •  Bessel{{ ),  b )  being  the  Bessel  distribution  of  mean 
0  and  variance  b.  The  surprisal  of  a(t)  is  by  definition: 

6{t)  =  u  (ot(t)) 

=  -In  (fa(a(t))) 

—  —  \n(A)  —  cos(a(£)),  V  t. 


(6.5) 
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Figure  6.2:  (a)  The  probability  density  function  of  a  von  Mises  distribution,  (b)  The 

corresponding  surprisal,  which  is  also  the  defined  manifold  J\f. 

6.4  Invariant  signature  for  planar  curves 

6.4.1  Definition 

In  [43],  Feldman  and  Singh  present  an  invariant  signature  for  planar  curves.  This 
invariant  is  based  on  the  turning  angle  <a(£),  the  change  in  the  direction  of  the  tangent  vector 
at  the  instant  t.  The  actually  considered  invariant  signature  is  the  information  gained  when 
measuring  a(t)  at  all  instants  t.  This  exactly  corresponds  to  0(t),  the  surprisal  of  a(t) 
as  defined  in  (6.5).  For  a  simple  planar  curve  of  length  L  sampled  at  N  equal  intervals 
of  arc  length  At,  a(t)  is  related  to  the  curvature  K,(t)  at  a  given  point  by  the  following 
approximation, 

a(t)  ~  K,(t)  •  At.  (6.6) 

We  note  that  scaling  the  curve  implies  scaling  both  n(t)  and  L  while  keeping  the  value  of 
tt(t)-^  invariant.  It  thus  follows  that  a(t)  is  a  measure  equivalent  to  /s(t)  except  that  it  is 
scale-invariant.  The  actually  considered  invariant  signature  is  the  information  gained  when 
measuring  a(t)  at  all  instants  £,  which  is  the  surprisal  function  9(t)  as  defined  in  (6.5). 

6.4.2  Illustration 

In  Figure  6.3  (a),  we  illustrate  an  example  of  a  planar  curve  (solid  blue)  and  its 
transformed  version  (dotted  red).  The  applied  transformation  is  a  combination  of  a  trans- 
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-1  —  Initial  planar  curve 


~*~  T ransformed  planar  curve _ _ 

-1.4  -1  -0.6  -0.2  0.2 

x 

(a) 


(b) 


Figure  6.3:  Two  similar  planar  curves  are  illustrated  in  (a).  In  (b),  the  same  turning  angle 
olt  is  obtained  for  the  two  planar  curves  in  (a). 


lation,  rotation,  and  scaling,  i.e.,  a  similarity  transform.  The  resulting  invariant  signatures 
0(t)  for  these  two  curves  are  exactly  the  same.  We  show  the  corresponding  turning  angle 
function  in  (Figure  6.3  (b)). 

6.5  Invariant  signature  for  space  curves 

6.5.1  Definition 

The  following  section  contains  the  major  contributions  of  this  paper.  We  start 
with  a  generalization  of  the  invariant  signature  we  presented  in  Section  6.4  from  2D  to  3D. 
In  our  case,  we  deal  with  space  curves  instead  of  planar  curves.  For  this  reason,  we  require 
two  turning  angles,  olt  and  ct#,  versus  one  in  the  planar  case;  hence,  we  naturally  use  these 
two  measurements  to  define  an  adequate  invariant  signature  curve,  olt  is  the  change  in  the 
direction  of  T.  Informally,  we  view  it  as  a  measure  of  how  much  a  curve  is  curved  and 
diverges  from  a  line,  olb  is  the  change  in  the  direction  of  B.  Similarly,  we  view  it  as  a 
measure  of  how  much  a  curve  is  twisted  and  diverges  from  a  plane. 

ar(t)  ~  n(t)  •  At  and  as(t )  ~  r{t)  •  At.  (6-7) 

Similarly  to  <a(£),  crr(^)  and  <a#(£)  follow  von  Mises  distributions  of  mean  zero  and  with  a 
spread  parameter  b  equal  to  1.  Thus,  exactly  as  in  (6.5),  these  two  measures  introduce  the 
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Figure  6.4:  From  the  surprisal  of  two  marginals  to  the  surprisal  of  one  binary  distribution. 


following  independent  invariant  signatures  for  space  curves: 


drit)  = 

and 

—  log(^4)  —  cos  (aT(t)) 

(6.8) 

Osit)  = 

-log(,4)  -  cos (aB(t)). 

(6.9) 

Instead  of  separately  considering  the  marginals  of  the  two  random  variables  ar(t)  and 
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Figure  6.5:  Manifold  defined  by  the  invariant  signature  curve  as  defined  in  Eq.(6.11). 


ob(£),  we  define  a  third  invariant  term  that  considers  the  random  vector  [cer(t), 

Thus,  the  distribution  fa  of  interest  becomes  the  binary  von  Mises  distribution  of  the 
independent  variables  ax  it)  and  such  that: 

=  A2  exp(cOS (o'j'(t))  +  cos (aB(t))).  (6.10) 

The  corresponding  surprisal  function  0(aT(t),  aB(t))  =  0(t)  becomes, 

0(t)  =  -In  (UipLTiOLB)) 

=  — 21n(^4)  —  cos(q't(^))  —  cos  (aB(t)),  (6.11) 

with  G  ([— 7r, tt])2. 
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Figure  6.6:  Two  different  sets  of  synthetic  space  curves.  In  (a),  T1  and  T2  are  two  similarity 
transforms.  In  (b),  (71,  (72  and  (73  are  three  similar  curves  except  at  one  inflection  point. 


6.5.2  Illustration  on  synthetic  space  curves 

We  test  the  invariant  signature  defined  in  (6.11)  on  some  synthetic  space  curves. 
In  Figure  6.6,  we  show  the  initial  curve  (71  defined  by, 

(71  :  t  E  [0 , 7r ]  — ►  (cos(2£), cos(3£), cos(5£)).  (6.12) 

We  define  two  similarity  transforms  through  the  matrices  T1  and  T 2;  therefore,  all  the 
curves,  except  (72,  are  similar  to  (71.  Indeed,  we  find  two  sets  of  invariants  as  shown  in 
Figure  6.7.  Those  corresponding  to  the  family  of  (71  are  represented  in  red  in  (a)  and 
(c)  (left  column),  (b)  and  (d)  correspond  to  (72  (right  column). 

This  example  illustrates  the  invariance  of  the  turning  angles  to  similarity  trans¬ 
forms.  It  is  actually  a  property  that  is  always  verified  thanks  to  the  intrinsic  nature  of 
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Turning  angle  of  tangent  vectors  T urning  angle  of  binormal  vectors 


(a) 

Turning  angle  of  tangent  vectors 


(c) 


(b) 

Turning  angle  of  binormal  vectors 


(d) 


Figure  6.7:  Invariants  for  the  space  curves  of  Figure  6.6.  (a)  and  (b)  correspond  to  turning 
angles  of  the  binormal  vectors,  (c)  and  (d)  correspond  to  turning  angles  of  the  tangent 
vectors.  In  red  are  the  signature  curves  for  C l’s  family.  In  green  are  the  ones  for  C 2. 
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Figure  6.8:  Two  different  sets  of  synthetic  space  curves.  In  (a),  T1  and  T2  are  two  similarity 
transforms.  In  (b),  C 1,  C 2  and  C 3  are  three  similar  curves  except  at  one  inflection  point. 


curvatures  (ft  and  r). 

6.5.3  Signed  angles  and  their  implementation 

We  test  a  tricky  case  where  local  versus  global  representations  are  confronted.  We 
use  the  new  curves  C 1,  C 2  and  C 3  illustrated  in  Figure  6.8,  and  defined  as  follows: 


Cl 

t  €  [— 7T,  7r]  - 

->  ( 4); 

(6.13) 

C2 

t  €  [— 7T,  7r]  - 

(i,sign(i)  x  t2, t4); 

(6.14) 

C3 

t  €  [ — 7T,  7r]  - 

->  (t,sign(i)  x  t2,sign(t)  x  t4); 

(6.15) 
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with 


sign(t)  = 


1  iff  t  ^  0; 
—  1  otherwise. 


Using  the  turning  angles  as  defined  in  Section  6.3,  we  find  the  turning  angles  represented 
in  Figure  6.9. 

Since  we  have  the  analytical  expression  of  the  curves  £71,(72  and  £73,  we  may 
compute  their  curvature  functions  ftl(t),  ^2(t),  and  ^3(t),  respectively,  as  well  as  their 
torsion  functions  rl(t),r2(t),  and  r3(t),  respectively.  We,  thus,  find 


and 


k1  (t) 

Tl(t) 


||(16t3, 12t2 ,  2)|| 

||(1,  2t,  4t3)  ||3  ’ 

48 1 

||(16t3,  — 12t2,  2) || 2  ’ 


(6.16) 

(6.17) 


k1  (t)  =  (t)  =  ^3(t); ;  (6.18) 

rl(t)  =  — sign(t)  x  r2{t)  =  sign(t)  x  r3(t).  (6.19) 


We  notice  that  these  curves  differ  in  their  torsion  functions.  Nevertheless,  the  turning 
angles  signatures  of  these  curves  fail  to  describe  their  difference.  This  phenomenon  of  local 
similarity  versus  global  dissimilarity  shows  the  necessity  of  taking  into  account  the  whole 
shape  representation.  This  is  equivalent  to  talking  about  a  memoryless  description  versus  a 
description  with  memory.  Therefore,  we  modify  the  initial  definition  of  the  turning  angles 
olt  and  olb •  The  idea  is  to  include  the  sign  negative  or  positive  for  our  turning  angles.  To 
that  end,  we  come  up  with  a  sign  convention  for  our  angles  in  space.  So  a  Frenet  frame 

at  a  time  t  is  now  considered  with  respect  to  the  Frenet  frame  at  time  ( t  —  St ),  or  in  our 

discretized  case  at  (t—  1).  Going  through  the  computations  and  simplifications  of  Appendix 
B,  we  present  the  following  definition  of  the  signed  turning  angles: 

aT(t)  =  sign(T(t)  •  N (t  -  1))  arccos(T(t)  •  T (t  -  1));  (6.20) 

aB(t)  =  sign(B(£)  •  T (t  —  1))  arccos(B(t)  •  B (t  —  1)).  (6.21) 

We  present  in  Figure  6.10  the  turning  angles  computed  using  (6.20)  and  (6.21).  In  Fig¬ 
ure  6.10  (a),  we  confirm  the  results  of  (6.18)  and  observe  a  perfect  correspondence  between 
the  turning  angles  corresponding  to  k  1,  /^2,  and  ^3.  We  recall  that  the  difference  between 
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(a) 


Figure  6.9:  Non  signed  turning  angles  for  the  space  curves  of  Figure  6.8. 
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discretized  arclength 


(a) 


(b) 


Figure  6.10:  Turning  angles  for  the  space  curves  of  Figure  4  (b). 


the  three  curves  is  in  the  sign  of  their  olt  functions  at  different  semi-open  intervals.  Using 
the  formulation  in  (6.21),  we  are  able  to  see,  experimentally,  in  Figure  6.10  (b),  overlaps 
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Figure  6.11:  Invariant  signature  curves  for  3D  modeling  curves  7I,  72  and  73;  the  same 
curves  appearing  in  Figure  6.5. 


and  symmetries  between  positive  and  negative  parts  of  the  turning  angles.  This  observa¬ 
tion  perfectly  complies  with  the  result  derived  in  (6.21).  This  change  in  the  sign  of  torsion 
functions  translates  the  effect  of  inflection  points  on  the  curves  Cl,  C 2,  and  C 3.  Actually, 
what  we  represented  in  Figure  6.10,  are  the  marginals  of  the  binary  distribution  presented 
in  (6.10).  The  actual  signature  curves  are  the  surprisal  of  this  distribution.  Their  accurate 
representation  should  be  on  the  curved  space  defined  in  (6.11).  We  show,  in  Figure  6.11, 
the  signature  curves,  for  Cl,  C2,  and  C3,  sitting  on  this  space.  In  Figure  6.10  (a)  and 
(b),  we  observe  overlaps  and  symmetries  between  some  parts  of  the  turning  angles.  This 
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observation  exactly  translates  what  is  happening  at  the  curves  level  because  of  the  effect  of 
the  inflection  points.  The  actual  signature  curves  sitting  on  a  subset  of  the  space  defined 
in  (6.11). 

6.6  Comparison  of  invariants 

The  choice  of  the  metric  to  use  on  any  model  is  very  important  in  assessing  the 
accuracy  of  this  one.  Indeed,  an  ineffective  metric  may  lead  to  a  failing  model  despite  the 
effectiveness  of  the  model  itself.  For  this  very  reason,  we  propose  to  walk  through  the  ideas 
that  led  to  the  final  metric  that  we  use  in  combination  with  the  proposed  representation 
for  space  curves. 

To  compare  two  given  curves  7 \  and  7 j,  we  compare  their  invariants.  In  this 
section,  we  explain  the  reason  behind  the  transition  from  (6.8)  to  (6.11).  Why  is  it  that  we 
used  the  binary  distribution  instead  of  the  marginals  even  though  the  two  variables  olt  and 
olb  are  assumed  to  be  independent  ? 

Let  us  assume  that  we  used  the  marginals  /a(cer)  and  In  such  case,  we 

approach  the  two  turning  angles  olt  and  olb  similarly  but  separately.  As  a  consequence, 
we  go  back  to  using  their  corresponding  surprisal  functions  6r(t)  and  9b  (t),  respectively, 
instead  of  9{t).  We  clearly  see  that  what  applies  to  0t  applies  to  9b]  hence,  we  proceed  by 
defining  a  metric  for  9t  and  expand  it  to  9b- 

If  the  curve  7 i  is  represented  by  N  sampling  points,  its  model  becomes  a  vector 
of  length  N.  Each  element  of  this  vector  ©j.(fc),  k  =  1,  •  •  •  ,  TV,  is  a  realization  of  9t  that 
takes  its  values  from  u([— 7r,7r]);  a  smooth  planar  curve  J\f  (or  a  1-dimensional  manifold)  as 
shown  in  Figure  6.14,  and  defined  as:  u  :  [—7 r,  tt]  — >  J\f  C  M2. 

To  compare  7^  and  7 j,  7 j  needs  to  be  sampled  with  exactly  the  same  number  of 
points  N.  That  is  7 j  would  be  assigned  a  similar  vector  0JT.  Comparing  7 i  and  7 j  in  a 
strictly  planar  case1  becomes  equivalent  to  comparing  the  two  vectors  QlT  and  0^. 

The  simplest  metric  to  use  would  be  a  mere  Euclidean  distance  between  the  two 
vectors  QlT  and  0^;  however,  it  is  easy  to  show  that  such  a  distance  would  be  inaccurate 
in  the  proposed  configuration.  Indeed,  if  we  look  at  the  example  of  Figure  6.12,  we  find 

1ln  3D,  in  addition  to  the  information  for  curvature,  one  needs  to  account  for  the  equivalent  information 
for  torsion,  z.e.,  some  vectors  ©h  and  0JB.  We  only  consider  the  planar  case  here  in  order  to  simplify  our 
explanations. 
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Figure  6.12:  Euclidean  versus  Riemannian  metric.  We  notice  that  for  two  turning  angles  oq 
and  0L2  that  are  the  same,  but  in  opposite  directions,  a  Euclidean  metric  cannot  distinguish 
between  their  corresponding  surprisals. 

that  for  two  turning  angles  aq  and  0L2  such  that  oq  —  —  <a2,  a  Euclidean  distance  clearly 
fails  to  distinguish  between  these  two  curves  as  it  would  give  a  null  distance.  It  is  therefore 
necessary  to  apply  a  Riemannian  metric  £(.,.)  on  the  manifold  J\f.  We  may  think  of  using 
a  normalized  geodesic  distance  between  two  realizations  Oi  and  Oj.  Specifically: 


(6.22) 


with  u(to)  =  Oi  and  u(t\)  —  Oj. 

While  this  metric  is  more  accurate,  we  still  have  to  ensure  that  the  two  curves 
7 i  and  7 j  are  sampled  with  exactly  the  same  number  of  points.  Enforcing  this  condition 
may  cause  correspondence  problems  as  shown  in  Figure  6.14.  One  may  think  to  deal  with 
this  problem  by  using  multiple  paths/comparisons  without  considering  or  trying  to  assign 
a  perfect  one-to-one  matching  between  sampling  points. 

Assuming  that  we  define  the  metric  that  does  the  work  for  the  comparison  of  ®lT  and  ©^, 
we  would  use  the  same  metric  to  compare  QlB  and  0JB.  A  new  question  would  still  arise; 
How  do  we  combine  the  two  metrics  £(©^,©^)  and  £(©#,©#)•  What  would  justify  a 
uniformly  weighted  average,  for  instance. 
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Figure  6.13:  Euclidean  versus  Riemannian  metric.  We  notice  that  for  two  turning  angles  oq 
and  0L2  that  are  the  same,  but  in  opposite  directions,  a  Euclidean  metric  cannot  distinguish 
between  their  corresponding  surprisals. 


In  order  to  avoid  any  arbitrary  considerations,  and  mainly  to  solve  all  correspon¬ 
dence  problems,  we  choose  to  use  the  binary  surprisal  function  6{t)  instead  of  the  marginals 
Oxit)  and  ds(t).  For  an  accurate  comparison  of  the  signature  curves,  we  define  a  Rieman¬ 
nian  metric  on  the  space  T  of  the  invariant.  The  invariant  we  defined  in  Eq.  (6.11)  is  a 
signature  curve  embedded  in  the  curved  space  T  created  by  the  two  variables  olt  and  olb 
(Figure  6.5).  All  the  invariant  signature  curves  that  we  are  to  compare  are  thus  constrained 
to  live  on  the  defined  invariant  space  that  we  call  T.  Defining  a  space  T  that  holds  all  the 
possible  invariants  is  a  natural  way  to  register  them,  and  therefore  deal  with  problems  of 
correspondence  between  signature  curves.  As  a  consequence,  we  may  directly  apply  a  dis¬ 
tance  measure  to  compare  these  invariant  curves  without  worrying  about  ensuring  a  prior 
registration.  We  thus  choose  to  compare  two  signature  curves  A1  and  A2,  corresponding  to 
two  space  curves  yl  and  y2,  by  considering  the  oriented  curve  Aa  =  A1  —  A2.  We  use  tools 
from  measure  theory  and  choose  to  refer  to  their  physical  intuition  in  relating  them  to  our 
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Figure  6.14:  Illustration  of  the  correspondence  problem 


problem  [45].  We  start  by  viewing  the  oriented  version  of  the  space  T  as  a  vector  field  F 
on  the  (27t  x  2n)  plane  defined  by  the  variables  olt  and  <ag.  We  directly  relate  F  to  V0(£), 
the  gradient  vector  field  of  0(£),  and  define  it  as  follows: 

F  :  ([— 7r,7r])2  -»•  M2 

(cxt,&b)  l_ ►  sin  (olt)  i  +  sin  (qlb)  j  • 

We  also  define  A*  ;  the  projection  of  Aa  on  the  (2n)2  plane.  A*  is  a  1-current  in  the 
space  dual  to  the  space  of  1-forms  P1([— tt,  n])2.  This  means  that  if  we  consider  F  (£)  = 
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F  (ar(t),  Q's(t))  and  for  all  0  from  D1([— tt,  7r])2, 


(6.23) 


(6.24) 


A 

With  these  notions  of  measure  theory,  we  naturally  use  the  flat  norm  F(y*  )  as  the  intrinsic 
distance  between  two  curves  yl  and  y2  whose  invariants  are  A1  and  A2,  respectively.  We 
thus  may  write, 


D 


where  A*  =  (a*1-A*2). 


F(A*1- A*2), 

sup{A*  (</>)  :  \\d(f)\\  ^  1  for  all  \\(j)\\  ^  1}, 


(6.25) 

(6.26) 


6.7  Correspondence  preservation 

In  what  follows  we  prove  and  simulate  different  properties  that  the  proposed  in¬ 
variants  along  with  their  metric  are  to  verify.  We  use  synthetic  space  curves  and  simulate 
different  scenarios. 

The  invariance  of  Euclidean  transforms  and  the  invariance  of  parametrization  are  naturally 
verified  from  the  properties  of  curvature  and  torsion.  We  also  showed  the  invariance  to 
scaling  in  Section  6.4.  We  therefore  start  by  looking  at  the  property  of  correspondence 
preservation.  Our  objective  is  to  show  the  advantages  of  registering  all  the  invariant  sig¬ 
nature  curves  on  the  space  T.  Some  of  these  advantages  are:  the  model’s  independence 
of  any  starting  point,  correspondence  preservation  in  occluded  curves  or  ensuring  a  partial 
matching  of  curves. 

Let  us  consider  two  space  curves  7I  and  72  represented  by  two  invariant  curves 
A1  and  A2,  and  defined  as  follows: 

71  :  [a,  b]  — >  M3;  (6.27) 

72  :  [c,d]^R3,  (6.28) 

such  that  a  ^  c  ^  d  ^  b  and  yl([c,  d])  =  /i(y2([c,  d])),  where  /i(-)  is  a  similarity  transform. 
In  other  words,  the  curve  y2  perfectly  matches  a  portion  of  7I  after  a  registration  operation, 
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(C) (d) 


(g) 


(h) 


Figure  6.15:  Illustration  of  the  problem  of  occlusion:  (a)  Initial  curve  7I.  (b)  Two  curves 
72  and  73  resulting  from  occluding  7I.  (c),  (d)  represent  the  turning  angles  for  7I.  (e),  (f) 
are  the  turning  angles  for  72,  and  (g),  (h)  for  73. 
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Figure  6.16:  Comparison  of  the  3D  objects  in  (a)  and  (c)  using  their  squigraph  represen¬ 
tations  in  (b)  and  (d),  respectively. 


i.e.,  /i_1(-),  in  R3. 

Since  the  invariant  curves  are  representations  that  are  locally  supported,  Aa  as  defined  in 
Eq.  (6.23),  becomes  the  invariant  curve  representing  yl ([a,  c])  and  7I ([&,  d])  that  we  call 
Al[a?c]  and  Al^p  The  distance  D ^71,72^  is  therefore  equivalent  to  F(Al[a?cj)  +F(A1^^). 
This  correspondence  problem  is  very  important  for  shapes  or  curves  that  are  subjected 
to  occlusion  or  provided  as  partial  views.  Let  us  for  instance  consider  the  curve  7I  in 
Figure  6.15  (a).  We  show  the  corresponding  turning  angles  olt  and  olb  in  Figure  6.15  (c) 
and  (d).  In  Figure  6.15  (b)  we  illustrate  the  same  curve  7I  but  this  time  subjected  to  an 
occlusion,  which  created  two  separate  curves  72  (in  red)  and  73  (in  green).  Computing  the 
turning  angles  for  each  one  of  these  two  curves  requires  using  two  distinct  starting  points; 
hence,  we  find  the  two  sets  of  signatures  (e)  and  (f)  for  72,  and  (g)  and  (h)  for  73.  It  is 
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Figure  6.17:  Invariant  for  partial  matching. 


clearly  difficult  to  directly  match  the  signatures  of  7I  with  those  corresponding  to  y3;  we 
face  a  correspondence  problem  caused  by  occluded  parts  of  a  curve. 

The  problem  becomes  less  challenging  when  using  the  invariant  space  T  as  a  way  to  register 
the  signatures.  Indeed,  Figure  6.16  presents  the  three  invariant  signature  curves  Al,  A2, 
and  A3  respectively  representing  7I,  y2  and  73. 

6.8  Partial  matching  of  primitive  shapes 

In  what  follows  we  give  an  example  of  an  partial  matching  of  3D  objects  (Fig¬ 
ure  6.18).  Using  the  modeling  technique  proposed  in  [10],  we  may  reduce  the  geometry  of 
a  primitive  3D  object  to  one  curve  in  R3  that  we  refer  to  as  a  modeling  curve;  thus  the 
full  vase  in  Figure  6.18  (a)  may  be  mapped  from  its  sampling  (Figure  6.18  (b))  to  just  one 
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curve  (Figure  6.18  (c)).  We  may  do  the  same  for  the  top  part  of  the  vase.  We  plot  the  two 
resulting  modeling  curves,  for  the  partial  and  the  full  vase,  in  Figure  6.18  (c).  We  notice 
that  it  is  nearly  impossible  to  see  any  partial  matching  by  simply  looking  at  these  curves. 
Instead,  we  extract  our  proposed  invariant  signature  curves.  We  find,  in  Figure  6.17,  a  nice 
partial  matching  of  these  signatures. 


6.9  Squigraph  comparison 


We  further  apply  these  signature  invariant  curves  to  compare  3D  shapes  through 
their  squigraph  representations.  The  particularity  of  these  skeletons  is  that  they  have  spatial 
curves  replacing  their  edges.  These  characteristic  curves  are  new  means  to  compactly  carry 
the  geometric  information  of  3D  shapes.  In  Figure  5.2  we  illustrate  a  typical  example  for 
which  the  practical  importance  of  the  proposed  invariant  signature  becomes  obvious.  The 
3D  shape  comparison  technique  simplifies  to  comparing  the  signatures  of  the  edge  curves 
with  the  same  colors  in  Figure  4.5.  For  these  skeletons,  we  define  a  new  global  metric  based 
on  Eq.  (6.25).  We  thus  consider  that  Cl  and  C 2  are  now  two  sets  of  curves  in  3D  such 
that  each  set  contains  N  curves  C\  and  C^,  i  =  1,  •  •  •  ,  iV,  respectively  representing  the 
geometrical  shapes  of  the  3-dimensional  parts  S\  and  i  =  1,  •  •  •  ,  iV,  that  constitute  each 
3D  object.  We  define  in  Eq.  (6.29)  the  new  distance  between  the  two  sets  Cl  and  C 2,  which 
is  also  the  distance  between  the  corresponding  objects  SI  and  S 2.  We  show  in  Chapter  7, 
the  results  of  using  this  distance  in  comparing  3D  objects. 


6.10  Conclusion 


In  this  Chapter,  we  presented  a  new  similarity  invariant  signature  for  space  curves. 
This  invariant,  since  based  on  the  tangent  and  the  binormal  turning  angles,  has  the  advan¬ 
tage  of  being  local,  unique  and  fully  invariant  to  similarity  transforms.  The  proposed 
invariant  proves  to  be  very  practical  to  use  in  3D  shape  modeling/comparison  problems. 
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Figure  6.18:  Comparison  of  the  3D  objects  in  (a)  and  (c)  using  their  squigraph  represen¬ 
tations  in  (b)  and  (d),  respectively. 
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Chapter  7 

Experimental  Results 


We  experimentally  investigate  the  performance  of  the  proposed  shape  modeling 
technique.  Throughout  our  definition  and  presentation  of  squigraphs,  we  have  shown  how 
rich  they  are  in  information.  In  what  follows,  we  illustrate  how  to  exploit  this  information 
for  different  levels  of  discrimination  and  for  different  application  objectives;  furthermore, 
we  investigate  the  robustness  properties  of  squigraphs  through  some  pertinent  examples. 

7.1  Discrimination  power 

To  assess  the  overall  discrimination  power  of  squigraphs,  we  compare  their  per¬ 
formance  to  those  of  well  established  approaches.  Our  comparison  involves  the  following 
techniques: 

•  Probability  density  function  (PDF)  descriptors:  Using  the  distributions  of  surface  fea¬ 
tures  is  a  technique  that  was  first  clearly  defined  and  analyzed  by  Funkhouser  et  al.  [50] . 
In  our  experiments,  we  use  the  GGF  as  the  surface  feature,  and  its  distribution  as  the 
shape  descriptor.  We  choose  the  Jensen- Shannon  Divergence  (JSD)  as  the  dissimi¬ 
larity  measure  between  the  GGF  distributions  [33,  48,  49,  50]. 

•  Classification  by  Characteristic  Resolution  (CCR)  [23,  51]:  We  have  shown  in  [23] 
that  each  class  of  shapes  determines  a  characteristic  resolution1.  This  parameter  is 
extracted  through  a  global  comparison  of  the  distribution  of  the  GGF  for  each  class  of 


1  Resolution:  is  the  number  of  vertices  used  to  describe  a  given  shape. 
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3D  objects.  In  short,  we  may  define  the  characteristic  resolution  as  being  the  lowest 
resolution  at  which  all  class  members  will  be  accurately  represented. 

•  Augmented  Multi-resolution  Reeb  Graph  (aMRG)  [2]:  This  technique  is  the  closest 
to  squigraphs  because  of  its  two-level  philosophy,  and  because  of  the  similar  Morse 
function  it  primarily  uses  in  defining  a  Reeb  graph.  aMRG  is  a  technique  that  has  been 
experimentally  established  to  outperform  many  3D  classification  techniques  including 
the  method  of  spherical  harmonics  [5]. 


Figure  7.1:  ROC  curves  comparing  different  classification  techniques. 

To  be  able  to  compare  all  these  different  techniques  on  the  same  basis  we  use  Receiver  Oper¬ 
ating  Characteristic  curves  (ROC)  as  illustrated  in  Figure  7.1.  We  then  use  the  Area  Under 
the  Curve  (AUC)  as  our  measure  for  classification  performance.  We  use  the  Princeton  [15] 
and  the  Technion  datasets  [52,  53,  54]  for  a  total  of  17  classes  and  239  objects.  The  overall 
performance  AUC  for  each  technique  is  summarized  in  Table  7.1.  With  a  little  more  than 
97%  of  overall  performance,  we  note  that  squigraphs  outperform  other  techniques.  Doing 
better  than  PDF  and  CCR  is  quite  expected  as  these  two  techniques  only  provide  a  global 
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Table  7.1:  Overall  performance  summary. 


Method 

PDF  CCR  aMRG  Squigraph 

AUC 

87.95%  92.53%  93.389%  97.35% 

description  of  a  complex  shape.  It  becomes  more  challenging  to  clearly  distinguish  between 
two  objects  whose  global  shapes  are  very  similar.  As  presented  earlier  in  Section  5.1,  aMRG 
and  squigraphs  are  multi-level  descriptors  which  enable  them  to  capture  global,  as  well  as 
local  features.  We  explain  the  better  performance  of  squigrpahs  over  aMRG  by  the  nature 
of  the  features  extracted  locally.  Indeed,  while  both  methods  extract  coarse  topological 
graphs  from  fully  intrinsic  representations,  i.e.,  a  geodesic  distance  independent  of  any  ref¬ 
erence  point,  aMRG  relies  on  local  features  that  relate  to  a  Euclidean-based  distance.  So, 
the  key  advantage  of  the  squigraph  representation  technique,  that  played  in  the  results  of 
Figure  7.1,  is  to  consistently  use  intrinsic  features  and  geodesic  distances. 

7.2  Primitive  shape  analysis 

In  order  to  thoroughly  understand  squigraphs  and  their  modeling  abilities,  we 
propose,  in  the  following  set  of  experiments,  to  investigate  the  different  properties  a  3D 
descriptor  is  required  to  enjoy.  Since  the  proposed  technique  proceeds  to  separately  consider 
homogeneous  parts  constituting  an  object,  i.e.,  mono-cardinality  subsurfaces,  it  is  intuitively 
easier  to  start  visualizing  the  different  properties  of  objects  with  primitive  geometries.  By 
primitive  geometry,  we  mean  that  the  corresponding  topological  graph  is  just  one  edge. 

7.2.1  Invariance  to  pose 

As  noted  earlier,  given  that  the  shape  of  an  object  remains  unaltered  when  sub¬ 
jected  to  a  similarity  transform,  we  require  a  geometric  descriptor  to  be  strictly  and  com¬ 
pletely  invariant  to  similarity  transforms.  To  demonstrate  the  invariance  of  a  GGF  to 
isometries  by  way  of  geometry  modeling,  we  carry  out  the  following  experiment.  Start¬ 
ing  with  the  shape  in  Figure  7.2  (a),  we  apply  each  of  the  following  transformations  to  it: 
(b)  rotation,  (c)  scaling,  (d)  shearing,  and  (e)  translation.  In  Figure  7.3  (a),  we  present  the 
resulting  five  modeling  curves  that  correspond  to  each  shape  of  Figure  7.2.  To  better  visu¬ 
alize  the  differences  of  the  curves  in  R3,  we  propose  to  project  them  on  the  plane  as  shown 
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Figure  7.2:  Invariance  to  pose.  The  original  object  represented  in  (a)  is  subjected  to: 
(b)  a  rotation  transform,  (c)  a  scaling  transform,  (d)  a  shearing  transform  and  in  (e)  to  a 
translation.  (Best  visualized  in  color) 


in  Figure  7.3  (b).  We  note  that  all  the  modeling  curves  show,  up  to  a  similarity  transform, 
a  similar  variation  to  that  of  the  curve  in  black  corresponding  to  the  initial  shape.  This 
experimental  result  is  consistent  with  the  results  of  Appendix  A. 

We  say  that  two  shapes  are  geometrically  identical  if  their  modeling  curves  are 
identical  up  to  a  similarity  transform.  For  this  reason,  we  use  the  representation  of  Chap¬ 
ter  6,  where  the  effects  of  similarity  transforms  are  eliminated  and  all  traces  corresponding 
to  the  modeling  curves  are  registered.  Identical  geometries  will  hence  see  their  traces  fully 
overlap.  We  present  the  traces  corresponding  to  this  experiment  in  Figure  7.4.  We  show 
the  top  view  of  these  traces  in  Figure  7.5.  These  traces  are,  by  definition,  invariant  to  all 
similarity  transforms.  We  thus  expect  to  see  the  result  of  Figure  7.5,  where  all  the  traces 
overlap  except  for  the  one  corresponding  to  the  shearing  transform;  moreover,  we  find  that 
the  distance  B  (as  defined  in  Eq.  (6.29))  between  the  original  shape  and  the  sheared  one 
is  equal  to  4.65  while  the  distance  of  other  shapes  away  from  the  original  one  is  between 
0.08  and  0.29.  Practically,  we  do  not  find  a  null  dissimilarity  between  geometrically  similar 
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Figure  7.3:  Modeling  curves  corresponding  to  the  objects  of  Figure  7.2.  (a)  Modeling 
curves  in  space,  (b)  Projection  of  the  modeling  curves  in  (a)  on  the  horizontal  plane. 


shapes  because  of  slight  differences  due  to  approximations  by  sampling  and  the  resulting 
Whitney  embeddings. 
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Figure  7.4:  Invariant  traces  for  curves  in  Figure  7.3. 


7.2.2  Progressive  deformation 


In  order  to  analyze  the  progression  of  the  distance  B,  as  a  result  of  progressively 
deforming  a  primitive  shape,  we  use  the  first  ellipsoid  starting  from  left  in  Figure  7.6  (a). 
This  ellipsoid  is  defined  by  the  following  equation: 

(x  -  X0 )2  (y  -  y0)2  ,  (z  -  z0 )2 


+ 


+ 


(rx)2  (ry)2  ( s  X  rz )2 

with  (x0,y0,z0)  =  (0,0,0),  (rx,ry,rz)  =  (|, \),  and  s  = 


=  1, 


(7.1) 


Each  time,  we  apply  a  vertical  shearing  transform  defined  by  the  matrix  T(s)  = 


/  1  0  0  \ 

0  1  0 

^  0  0  S  j 

s  being  the  shearing  factor.  We  vary  s  from  |  to  3.  We  denote  the  resulting  ellipsoid  for 
a  given  s  by  £(s).  We  note  that  £(1)  is  the  2-dimensional  sphere  centered  at  the  origin  (0, 


76 


Figure  7.5:  Invariant  traces  for  curves  in  Figure  7.3. 


0,  0)  with  a  radius  equal  to  Before  proceeding  with  the  analysis  of  the  present  experi¬ 
ment,  we  have  to  answer  a  very  important  question  that  arises  when  considering  singular 
shapes  such  as  spheres.  We  first  recall  that  we  only  talk  about  a  Morse  function  when  we 
ensure  that  all  critical  points  are  non-degenerate;  thus,  the  GGF  on  a  spherical  shape  is 
not  a  Morse  function.  Our  proposed  solution  is  to  slightly  disturb  such  a  singular  shape. 
For  a  spherical  shape  for  instance,  we  may  disconnect  a  polar  point  from  other  points.  As 
illustrated  in  Figure  7.7,  the  distribution  of  the  GGF  on  a  spherical  surface  quickly  changes 
from  a  simple  dirac  function  to  the  distribution  shown  in  Figure  7.7  (d). 

We  apply  this  same  perturbation  technique  to  all  the  shapes  in  Figure  7.6  (a),  and  extract 
a  modeling  curve  for  each  set  of  iso-geodesic  curves  in  Figure  7.6  (b).  Figure  7.8  illustrates 
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Figure  7.6:  Comparison  of  the  geometry  of  five  deformed  spheres,  (a)  shows  the  5  consid¬ 
ered  shapes,  and  (b)  shows  the  correspondingh  spatial  modeling  curves. 


the  resulting  modeling  curves.  Our  first  observation  is  that  these  curves  are  planar,  i.e.,  in 
R2,  while  the  proposed  Whitney  modeling  technique  is  in  R3.  This  is  due  to  the  simplicity 
of  the  evolution  of  the  iso-geodesic  curves.  Very  briefly,  we  state  that  more  variations  imply 
more  dimensions. 

Focussing  on  this  important  aspect  of  predicting  and  understanding  the  modeling  curves 
constitutes  a  future  research  direction.  Indeed,  the  present  experiment  and  its  results  lead 
to  further  investigations  toward  defining  the  applicability  and  limitations  of  the  strong 
Whitney  embedding  versus  the  easy  Whitney  embedding  theorem2. 

Our  choice  of  a  simple  progressive  vertical  shearing  transform  is  motivated  by  two  points: 
first,  applying  a  simple  directional  geometric  deformation,  and  second,  being  able  to  quan¬ 
titatively  describe  this  deformation.  In  the  present  example  we  use  the  shearing  factor  s  for 
this  quantification;  hence,  we  are  able  to  visualize  the  variation  of  B  versus  s;  thus,  the  red 

2 Indeed,  while  the  easy  Whitney  embedding  theorem  allows  embedding  an  n-dimensional  Hausdorff  man¬ 

ifold  into  R2n+1,  the  strong  Whitney  embedding  allows  going  lower  and  defining  an  embedding  in  R2n  [38]. 
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Figure  7.7:  The  original  sphere  represented  in  (a)  is  disturbed  to  create  a  non-degenerate 
GGF  function,  whose  color  mapping  is  shown  in  (b).  The  corresponding  distributions  of 
the  GGF:  (c)  Before,  and  (d)  after  disturbing  the  sphere. 
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Figure  7.8:  Analysis  of  the  effect  of  shearing  transforms  on  the  modeling  curves. 


full  curve  in  Figure  7.9  shows  the  smooth  evolution  of  the  distance  between  £(\)  and  other 
ellipsoids  with  |  <  s  ^  3.  Taking  a  sphere  as  the  reference,  i.e.,  £(1),  we  may  see  the  effect 
of  the  vertical  deformation  in  the  two  directions  s  <  1  and  s  >  1.  The  resulting  distance 
curve  is  illustrated  in  Figure  7.9  with  a  dotted  blue  line.  We  note  that  the  distance  B  follows 
an  exponential  behavior  as  it  starts  slowing  down  for  large  values  of  As.  One  may  explain 
this  tendency  by  referring  back  to  the  human  perception.  In  fact,  looking  at  Figure  7.6  (a), 
we  may  classify  these  ellipsoids  into  two  categories:  category  1,  corresponding  to  s  <  1;  and 
category  2,  corresponding  to  s  >  1,  the  sphere  (s  =  1)  representing  the  transition  point.  We 
then  may  say  that  when  an  ellipsoid  is  deformed  to  be  away  from  its  category,  the  distance 
B  is  relatively  high,  but  once  it  reaches  a  new  category,  additional  deformations  in  the  same 
direction  will  only  add  slight  distances.  We  find,  for  instance,  that  B(£(2),£(3))  =  0.045 
while  B(£ (0.5),  £(1.5))  =  0.22,  that  is  about  5  times  the  first  distance  and  for  the  same  As. 
The  present  result  constitutes  an  important  step  in  understanding  the  geometry  of  shapes, 
and  translating  the  human  perception  of  geometry. 
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Figure  7.9:  Analysis  of  the  effect  of  shearing  transforms  on  the  modeling  curves. 


7.2.3  Mono-cardinality  subsurfaces 

In  this  experiment  we  continue  our  observation  and  analysis  of  mono-cardinality 
subsurfaces,  but  this  time  we  consider  more  complex  shapes  found  in  real-world  objects. 
We  thus  compare  selected  parts  that  constitute  the  edges  of  the  squigraphs  extracted  for 
the  datasets  used  in  Section  7.1. 

We  summarize  our  comparison  results  with  the  confusion  matrix  in  Figure  7.10.  We  note 
that  when  legs  are  bent,  the  modeling  curves  are  still  able  to  detect  a  difference  between 
stretched  and  bent  legs.  This  phenomenon  is  equivalent  to  the  example  of  Subsection  6.5.3 
where  we  illustrated  how  the  modeling  curves,  and  squigraphs  in  general,  are  able  to  detect 
an  inflection  point.  Depending  on  the  application,  one  might  think  of  bent  and  stretched 
legs  (or  arms)  as  being  the  same.  In  this  case,  we  note  that,  while  modeling  curves  are  still 
able  to  be  lenient  on  the  effect  of  one  inflection  point,  we  recommend  for  simplicity  and 
efficiency  to  directly  use  the  distribution  of  the  GGF  along  each  part. 
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Figure  7.10:  Matching  and  comparison  of  the  geometry  of  mono-cardinality  subsurfaces. 
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Figure  7.11:  Comparison  and  matching  of  6  different  bipedal  subjects. 


7.2.4  Shape  matching 

The  synthesis  of  the  distances  between  different  mono-cardinality  subsurfaces  pro¬ 
vides  a  global  distance  for  complex  shapes  as  given  by  Eq.  6.29.  We  present  a  couple  of 
relevant  examples  of  confusion  matrices  for  different  classes  and  subclasses  of  3D  shapes.  In 
Figure  7.11,  we  note  the  power  of  the  squigraph  technique  in  differentiating  between  objects 
that  are  topologically  similar,  but  who  belong  to  distinct  object  classes.  In  Figure  7.12, 
we  present  the  result  of  comparing  the  same  object /subject,  but  in  different  poses.  We 
note  again  the  ability  of  squigraphs  in  finding  differences  due  to  both  rigid  and  non-rigid 
transformations.  The  best  way  to  understand  these  results  is  to  go  back  to  the  analysis  of 
changes  that  occur  on  primitive  shapes.  In  this  case,  it  means  observing  what  happens  to 
each  part  of  a  human  body  separately. 
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Figure  7.12:  Comparison  and  matching  of  the  14  different  poses  for  the  same  subject 
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Figure  7.13:  Robustness  to  noise  by  short  secants  filtering. 


7.3  Robustness  evaluation 

For  the  practical  applicability  of  the  modeling  space  curves  and  the  final  squigraphs, 
it  is  important  to  evaluate  their  robustness  to  both  noise  and  object  discrete  decimation. 

7.3.1  Robustness  to  noise 

As  for  any  practical  problem,  the  measurements,  when  scanning  3D  models,  con¬ 
tain  additive  noise.  The  secant-based  method  (described  earlier),  being  directly  dependent 
on  the  data,  raises  natural  concern  about  the  level  of  accuracy  of  the  modeling  curves  in 
presence  of  noisy  measurements.  For  testing  purposes,  we  add  Gaussian  noise  with  an  am- 
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Figure  7.14:  Secants  filtering  with  different  thresholds 


plitude  ranging  between  0  to  1%  of  the  shape  bounding  radius  for  each  considered  object. 
Figures  7.13  (a)  and  (b)  illustrate  the  object  “Vase”  with  and  without  additive  noise.  In  Fig¬ 
ure  7.13  (e),  we  observe  the  impact  of  noise  on  the  modeling  curve.  The  dissimilarity  score 
between  the  initial  shape  “Vase”  and  its  noisy  version  is  of  3.08.  This  basically  indicates 
that  the  two  modeling  curves  are  very  different  and  may  induce  an  erroneous  recognition 
result.  Broomhead  and  Kirby  [40]  propose  a  solution  to  alleviate  the  effect  of  noise  on 
embeddings.  They  show  that  the  shortest  secants  induce  the  most  severe  distortion.  We 
propose,  as  a  result,  to  filter  the  shortest  secants  prior  to  applying  the  LTMADS  algorithm. 
We  test  the  effect  of  filtering  out  the  short  secants  where  the  threshold  is  a  fraction  of  the 
longest  secant.  Figure  7.14  illustrates  the  result  of  250  Monte  Carlo  simulations.  It  is  clear 
that  the  errors  due  to  noise  are  reduced  after  filtering,  with  a  more  dramatic  impact  for 
larger  thresholds. 
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Figure  7.15:  Robustness  to  decimation. 


7.3.2  Robustness  to  decimation 

The  resolution  of  a  3D  object  measurement /scanning  may  vary  with  applications, 
and  it  is  hence  important  to  evaluate  our  proposed  geometric  modeling  in  presence  of 
such  variations.  To  that  end,  we  consider  an  objet  “Frog”  to  assess  the  robustness  of  the 
modeling  to  decimation.  Taking  the  highest  resolution,  of  10000  vertices,  as  a  reference,  we 
measure  the  dissimilarity  induced  by  a  progressive  decrease  in  resolution  for  object  “Frog”. 
We  illustrate  this  experiment  in  Figure  7.15.  As  explained  in  [23],  using  the  GGF,  we 
may  assign  a  particular  resolution,  to  each  shape,  referred  to  as  a  characteristic  resolution. 
By  definition,  the  characteristic  resolution  corresponds  to  the  lowest  number  of  vertices 
able  to  accurately  represent  a  considered  shape.  For  the  shape  “Frog”,  we  find  that  the 
characteristic  resolution  is  about  5000  points.  Analyzing  the  results  of  Figure  7.15,  we  may 
conclude  that  for  “Frog”,  our  modeling  is  accurate  for  resolutions  greater  or  equal  to  the 


87 


characteristic  resolution  of  5000  points.  Thus,  our  geometric  modeling  is  consistent  and 
nearly  invariant  to  decimation. 

7.4  Conclusion 

In  this  Chapter,  we  presented  a  novel  and  effective  geometric  modeling  technique 
which  significantly  reduces  a  3D  shape  representation  to  a  squigraph.  The  first  positive  re¬ 
sults  are  encouraging  to  warrant  further  study  and  more  expanded  applications.  The  overall 
objective  of  our  work  is  to  develop  a  comprehensive  and  efficient  classification  algorithm 
which  can  further  translate  the  human  understanding  of  the  geometry  of  shapes  from  a 
local  to  a  global  viewpoint. 
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Chapter  8 

Adaptive  Embedding  and  its 
Application  to  Network  Failure 
Detection 


We  define  a  new  adaptive  embedding  approach  for  data  dimensionality  reduction 
applications.  Our  technique  entails  a  local  learning  of  the  manifold  of  the  initial  data, 
with  the  objective  of  defining  local  distance  metrics  that  take  into  account  the  different 
correlations  between  the  data  points.  We  choose  to  illustrate  the  properties  of  our  work  on 
the  isomap  algorithm.  We  show  through  multiple  simulations  that  the  new  adaptive  version 
of  isomap  is  more  robust  to  noise  than  the  original  non-adaptive  version;  moreover,  we  use 
our  proposed  technique  to  detect  intrusions  in  sensor  networks. 

8.1  Introduction 

In  recent  years,  data  sizes  have  drastically  increased.  As  a  result,  there  has  been 
a  great  research  focus  on  improving  and  defining  effective  dimension  reduction  techniques. 
These  efforts  are  extremely  relevant  if  not  crucial  to  data  storage,  visualization,  and  analysis 
application.  The  growing  demand  of  storage  and  archiving  resources,  together  with  their 
inefficient  current  exploitation,  and  their  increasing  cost,  have  also  made  this  line  of  research 
very  relevant.  The  objective  behind  learning  and  reducing  the  dimension  of  data  is  to 
eliminate  any  redundant  information  while  still  preserving  the  intrinsic  and  underlying 
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structure.  One  may  think  of  this  problem  as  an  attempt  to  find  all  the  variables  that 
may  be  combined  into  fewer  variables  without  destroying  the  interactions  between  the  data 
points.  To  formulate  this,  we  proceed  as  follows. 

8.1.1  Problem  statement 

Given  a  data  point  sample  of  N  points  i  =  1,  •  •  •  ,  N,  from  a  d-dimensional 
smooth  manifold  J\4,  where  M  C  and  d  <  n,  in  a  dimension  reduction  problem,  we  aim 
at  finding  a  homeomorphism  /(  • ),  or  the  image  of  {Xi\  by  a  homeomorphism  /(  •  )  such 
that: 

f  :M  Cl"  -*■  Rd  (8.1) 

Xi  »  f(Xi)=Yi,  for  i  =  l,-.-,N. 

The  goal  is  to  effectively  determine  the  data  {Yi}  of  reduced  dimension  d  that  reflects 
all  the  information  present  in  the  initial  data  {Xi}.  Explicitly  finding  the  function  /(•) 
is  however,  not  required  in  dimension  reduction  applications,  as  the  approach  is  in  spirit 
non-parametric  and  data  driven.  We  may  distinguish  two  classes  of  dimension  reduction 
approaches;  The  first  class  includes  all  the  classical  methods,  or  linear  methods,  such  as 
Principal  Component  Analysis  (PCA)  [55],  and  Multidimensional  Scaling  (MDS)  [56].  In 
contrast,  the  second  class  corresponds  to  non-linear  techniques  [58],  also  referred  to  as  man¬ 
ifold  learning  methods.  The  linear  algorithms  are  well  known  and  fully  understood.  They 
are  often  used  to  preprocess  data  before  proceeding  to  further  analysis.  The  fundamental 
idea  of  PCA,  a  widely  used  linear  algorithm,  is  to  only  keep  predominant  dimensions  of 
the  data.  Geometrically  speaking,  this  corresponds  to  linearly  projecting  the  data  on  the 
vectors  of  highest  variances;  however,  this  only  works  if  the  data  has  a  linear,  or  close  to 
linear,  structure.  Any  nonlinearity  in  the  data  can  only  be  accounted  for  by  more  general 
methods  extending  beyond  linear  structure.  The  second  class  of  manifold-based  algorithms 
tends  to  find  the  same  predominant  feature  axes  using  more  relaxed  assumptions.  The 
common  assumption  of  these  more  powerful  techniques,  is  that  all  the  data  samples  live  on 
a  Riemannian  manifold  embedded  in  a  high  dimensional  Euclidean  space.  These  techniques 
have  proven  so  far,  to  be  the  most  successful  approaches  [57].  When  the  source  of  data 
has  a  relatively  few  degrees  of  freedom,  it  becomes  possible  to  reduce  the  dimensionality  of 
this  data  while  keeping  its  inherent  structure.  In  fact,  a  lot  of  real  world  applications  fall 


90 


within  this  class  of  underlying  lower  dimensional  systems  which  may  easily  be  formulated 
in  a  manifold  learning  setting.  Such  applications  include  face  recognition  [68] ,  face  pose  de¬ 
tection  [65],  gait  analysis  [66],  and  human  motion  data  interpolation  [67].  There  are  about 
four  widely  known  manifold  learning  techniques;  Locally  Linear  Embedding  (LLE)  [59], 
Laplacian  eigenmap  [60,  61],  Hessian  eigenmap  [62],  and  isomap  [63].  Manifold  learning 
algorithms  always  assume  the  observed  cloud  of  data  points  as  part  of  a  smooth  manifold. 
To  proceed  with  the  analysis  of  this  data,  we  start  by  constructing  a  graph  connecting  all 
the  data  points  and  preserving  the  structure  of  the  manifold.  One  usually  defines  an  e  —  ball 
neighborhood  of  fixed  radius  around  each  data-point  to  carry  out  the  analysis.  All  these 
techniques  have  shown  very  successful  results  in  ideal  conditions;  nevertheless,  there  is  very 
limited  work  in  addressing  the  effect  of  noise  and  the  choice  of  the  neighborhood  size.  Very 
simple  experiments  may  show  how  crucial  it  is  to  take  these  considerations  into  account. 

Our  goal  in  this  work  is  to  address  the  noise  problem,  and  propose  a  way  to  develop 
a  new  manifold-learning  technique,  with  a  built  in  robustness  to  noise.  The  key  idea  is  to 
replace  the  arbitrary  choice  of  an  approximate  Euclidean  distance,  and  to  instead  use  a 
locally  adaptive  distance.  To  achieve  that,  we  propose  to  account  for  sample  data  points’ 
correlations  in  defining  their  neighborhood.  We  follow  the  approach  of  Carlsson  et.  al  in 
addressing  the  problem  of  an  “adequate”  neighborhood  size  [72].  So,  to  go  beyond  the  one 
neighborhood  size  which  achieves  some  particular  embedding,  and  to  gain  further  insight, 
we  argue  that,  geometrically  and  topologically,  there  is  more  interest  to  discover  by  evolving 
a  neighborhood  size.  We  demonstrate  in  addition  the  flexibility  of  this  approach  in  a  sensor 
network  intrusion  detection  application. 

The  remainder  of  the  paper  is  organized  as  follows:  Section  8.1.2  is  an  overview  of  the 
different  existing  approaches  to  robust  manifold  embedding  techniques,  with  more  emphasis 
on  the  isomap  algorithm.  In  Section  8.2,  we  discuss  the  classical  isomap  algorithm  which,  in 
contrast  to  our  proposed  technique,  is  non-adaptive.  In  Section  8.3,  we  describe  in  detail  our 
proposed  adaptive  method.  We  evaluate  the  benefit  of  an  adaptive/non- adaptive  isomap 
in  Section  8.4,  using  a  residual  covariance  as  a  performance  measure.  We  illustrate  our 
proposed  technique  by  analyzing  synthetic  data  in  Section  8.5,  and  demonstrate  its  wider 
applicability  by  solving  a  sensor  network  failure  detection. 
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8.1.2  Related  work 

Many  of  the  existing  manifold  learning  techniques  show  successful  results  on  some 
well  chosen  data;  and  they  all  share  a  limiting  failure  when  in  presence  of  more  challenging 
datasets.  The  difficulties  are  often  due  to  the  intrinsic  topological  and  geometric  structure 
of  the  manifold  with  quick  variations  in  their  curvature  and  non-convex  boundaries  [69]. 
Additional  difficulties  result  from  the  properties  of  the  real  data  such  as  the  sampling  dis¬ 
tribution  and  the  nature  and  level  of  the  prevailing  noise. 

To  address  existing  limitations  and  to  further  improve  the  embedding  results,  and  extend 
the  applicability  of  the  current  manifold  learning  techniques,  we  propose  to  account  for 
these  overlooked  characteristics.  The  idea  is  to  progressively  adapt  to  the  data  at  hand 
in  tracing  the  local  connectivity  between  the  point  samples.  Some  recent  efforts  have  ex¬ 
plored  adaptive  manifold  learning,  by  specifically  focusing  on  two  parameters:  the  intrinsic 
dimensionality  of  the  data,  and  the  size  of  the  neighborhood.  Wang  et  al.  propose  in  [69]  a 
method  to  adaptively  select  the  neighborhood  size.  They  base  their  technique  on  determin¬ 
ing  the  alignment  space  of  local  tangents.  Cost  et  al.  on  the  other  hand  define  an  intrinsic 
dimensionality  using  K- nearest  neighbors  graphs  [70].  In  [71],  Levina  et  al.  adopt  a  local 
estimate  for  the  intrinsic  dimension  at  each  point. 

Our  present  work  also  focusses  on  an  adaptive  embedding  of  the  manifold;  we,  in  contrast, 
point  out  an  additional  characteristic  that  appears  to  be  as  critical  as  the  choice  of  the 
embedding  dimension  or  the  neighborhood  size.  We  indeed  show  in  what  follows,  that  the 
choice  of  a  Euclidean  distance  is  suboptimal  in  determining  the  local  connectivity  between 
data  points,  and  therefore,  introduce  a  new  adaptive  distance  locally  defined  for  each  point. 

8.2  Non-linear  manifold  searching  techniques 

As  noted  earlier,  our  proposed  effort  builds  on  existing  manifold  embedding  tech¬ 
niques;  therefore,  we  start  by  recalling  the  preliminary  steps  of  a  non-linear  manifold  learn¬ 
ing  technique. 

8.2.1  Isometric  feature  mapping  (Isomap) 

Among  the  multiple  manifold  learning  algorithms,  we  choose  to  use  the  isomap 
algorithm  to  illustrate  our  ideas.  This  choice  is  due  to  the  isomap  success  in  numerous 
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Figure  8.1:  Different  phases  of  isomap. 


embedding  problems  and  its  well  established  properties  [63].  The  principle  and  motivations 
of  this  work  are,  nevertheless,  extendable  to  other  mapping  (embedding)  algorithms. 
Isomap  is  a  nonlinear  mapping  algorithm  which  starts  with  the  assumption  that  the  high 
dimensional  data  lie  on  a  Riemannian  manifold  [64].  To  achieve  the  dimension  reduction, 
isomap  defines  a  mapping  that  aims  to  preserve  the  geodesic  distances  on  the  initial  mani¬ 
fold.  We  may  describe  isomap  as  merely  an  improved  version  of  Multidimensional  Scaling 
(MDS)  embedding  where  the  interpoint  distance  is  a  geodesic,  i.e.,  restricted  to  lie  on  the 
initial  manifold  of  the  data. 

Concretely,  for  a  data  point  sample  of  N  points  X^i  —  1,  •  •  •  ,  TV,  from  a  d-dimensional  man¬ 
ifold  M,  where  M  C  and  d  <  n,  we  detail  the  different  steps  of  the  isomap  embedding 
algorithm  in  Table  8.1. 

8.2.2  Local  connectivity  graph 

The  very  key  idea  of  isomap  is  to  consider  geodesic  distances  along  the  manifold 
of  data.  To  practically  approximate  the  intrinsic  geodesic  distance  on  a  manifold,  we  need 
to  locally  connect  each  point  to  its  k  nearest  neighbors,  or  equivalently  to  the  neighbors 
within  the  e- neighbor  hood.  By  so  doing,  we  result  in  a  graph  that  approximates  the  real 
manifold.  In  Figure  8.1  we  illustrate  the  different  phases  of  isomap;  from  a  continuous 
manifold  Figure  8.1  (a)  to  its  approximating  graph  Figure  8.2  (a),  and  finally  the  embedding 
Figure  8.2  (b).  We  show  in  Figure  8.2  (b)  how  severely  a  connectivity  graph  may  be  affected 
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Table  8.1:  Non-adaptive  isomap  algorithm. 

Step  1  :  Construct  a  weighted  graph  G; 

Let  G  =  {A,  X},  where:  m  A  is  an  (N  x  N)  adjacency  matrix; 

>X=[X1,---,XJV]T; 

Compute  D e,  the  matrix  of  Euclidean  distances  between  each  two  points  in  {Xj}£Li- 
Choose  e,  the  neighborhood  radius. 

for  i,  j  G  {1,  •  •  •  ,  N}  do 

if  D#(i,  j)  <  e  do 

A(i,j)  =  D  E(,i,j); 

else 

A(i,j)  =  oo; 

end  if 

end  for 

Step  2:  Compute  geodesic  distances  on  G. 

Let  T>g  be  the  matrix  of  geodesic  distances  between  each  two  points  in  {Xi\f=1. 
do  Dg  =  A;  (initialization) 
for  z,  j  G  {1,  •  •  •  ,  iV };  k  =  1;  do 

while  D G(i,j)  ^  DG(i,k)  +  T>G(kJ)  do 
for  k  G  {1,  •  •  •  ,  TV}  do 

Dg(l  j)  =  min(DG(i,  j),DG(i,  fc)  +  DG(fe,  j)); 

end  for 
end  while 

end  for 

Step  3:  Apply  MDS  on  DG. 


in  the  presence  of  noise.  This  consequently  yields  an  inaccurate  embedding  of  a  given 
manifold  as  shown  in  Figure  8.2  (d). 

8.3  New  adaptive  distance 

The  graphs  illustrated  in  Figure  8.2  (b)  and  (d)  are  the  result  of  considering  a 
Euclidean  neighborhood.  We  herein  maintain  that  the  choice  of  the  distance  is  crucial 
in  constructing  good  connectivity  graphs.  Our  objective  is  to  define  a  more  appropriate 
distance  that  alleviates  the  effect  of  noise  and  obtain  accurate  graphical  approximations. 
In  what  follows  we  provide  the  intuitive  rationale  for  the  choice  of  a  new  adaptive  distance. 
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Figure  8.2:  Failure  of  isomap  in  a  noisy  setting. 


We  subsequently  present  a  mathematical  formulation  of  new  solutions  to  the  embedding 
problem  to  result  in  an  improved  technique  described  in  Table  8.1. 
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Figure  8.3:  Mahalanobis  vs  Euclidean  distance:  On  the  left,  the  result  of  using  a  Euclidean 
distance.  On  the  right,  the  result  of  using  a  Mahalanobis  distance. 


8.3.1  Motivation 

In  spite  of  the  isomap  good  embedding  results,  it  remains  very  unstable  and  sensi¬ 
tive  to  noise,  as  well  as  to  the  choice  of  the  parameter  e  and  the  distance  function  used  prior 
to  applying  MDS.  Changing  the  distance  from  Euclidean  to  geodesic  thus  appears  to  be  in¬ 
sufficient  to  completely  preserve  the  intrinsic  geometric  structure  of  the  initial  manifold  M. 
In  what  follows,  we  propose  to  account  for  the  statistical  properties  of  the  observed  data. 
Specifically,  our  technique  consists  in  considering  the  correlation  between  each  point  and 
the  rest  of  the  observed  data  points,  and  subsequently  exploit  this  information  to  connect 
it  to  its  neighbors.  This  idea  is  exactly  equivalent  to  using  a  Mahalanobis  distance  [73].  To 
better  understand  the  intuition  behind  our  choice,  we  illustrate  the  result  of  constructing 
a  graph  connectivity  for  the  sample  points  in  Figure  8.3.  We  note  that  when  using  a  Eu¬ 
clidean  distance  to  determine  the  neighbors  of  each  data  point  causes  a  miss  of  some  details 
in  the  structure  of  the  data  set. 

8.3.2  Adaptive  embedding  algorithm 

In  what  follows  we  propose  a  new  distance  matrix  Dm  to  replace  D#  in  the 
algorithm  described  in  Table  8.1.  Our  objective  is  to  define,  each  time,  a  distance  that  is 
fully  dependent  on  the  sample  point  hence,  we  rescale  the  data  coordinates  based 
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Table  8.2:  Description  of  the  learning  step  (  new  adaptive  distance) 
Step  0:  Compute  Dm,  the  new  distance  on  M. 


Choose  e i,  the  neighborhood  radius  for  manifold  learning; 
and  62,  the  neighborhood  radius  for  the  construction  of  G. 

for  i  G  {1,  •  •  •  ,  N}  do 

=  Xj\  (initialization) 
for  j  =  1,  cdots ,  N  do 

while  D#(i,  j)  <  ei  do 
Y*  =  [Yjj  Xj']; 
end  while 

end  for 

JT  —  cov(Y^),  cov(-)  being  the  covariance  matrix; 

end  for 

for  i  G  {1,  •  •  •  ,  N}  do 

for  j  G  {1,  •  •  •  ,  N}  do 

D  =  (V-  -  V)  -  Xif; 

end  for 
end  for 

do  e  =  e2;  =  Dm; 

go  to  Step  1.  (See  Table  8.1) 


on  their  distributions  on  M  as  well  as  their  correlations.  Since  this  technique  relies  on 
a  learning  procedure  and  directly  uses  isomap  to  build  on,  we  refer  to  it  as  an  adaptive 
isomap  algorithm.  We  hence  use  the  algorithm  of  Table  8.1  with  a  learning  step,  i.e,  Step 
0,  as  described  in  Table  8.2. 

8.4  Performance  comparison 

In  this  section,  we  qualitatively  and  quantitatively  compare  the  performances  of 
the  two  versions  (adaptive  and  non-adaptive)  of  isomap  embeddings.  To  that  end,  we  start 
by  defining  a  performance  measure.  We  subsequently  simulate  different  classical  examples 
of  manifolds  to  embed  in  a  lower  dimensional  space.  The  choice  of  our  examples  is  such 
that  one  may  visually  inspect  and  verify  the  properties  as  well  as  the  intuition  behind  each 
technique. 
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8.4.1  Residual  variance 

We  choose  a  residual  variance  p  to  be  our  performance  indicator.  In  further  ap¬ 
plications,  we  may  use  it  to  investigate  the  topological  structure  of  a  manifold.  Residual 
variance  p  is  defined  in  Eq.  (8.2). 

p—  1  —  (corrcoef(Dx,  Dz))2,  (8.2) 


where: 

•  Dx:  is  the  distance  matrix  for  the  initial  data  in  M.  For  an  isomap  embedding 
technique,  this  matrix  is  the  geodesic  distance  matrix  which  is  exploited  by  MDS. 
We  note  that  may  vary  as  it  depends  on  the  first  distance  matrix  used  to  identify 
a  neighborhood.  For  the  classical  non-adaptive  isomap,  this  distance  was  simply 
Euclidean,  i.e.,  T>e,  while  for  the  adaptive  case  we  define  it  as  D m-  This  difference  is 
crucial  to  comparing  the  performance  improvement  of  the  adaptive  isomap  algorithm 
over  its  classical  counterpart. 

•  Dz:  is  the  distance  matrix  for  the  final  (embedded)  data  of  reduced  dimension  p 
(p  <  n).  We  take  advantage  of  the  simplicity  of  the  geometries  of  our  examples  (swiss 
rolls,  hemispheres,  parallel  sheets),  and  take  p  equal  to  2  and  consider  Dz  to  be 
Euclidean.  It  becomes  trivial  to  visually  verify  the  accuracy  of  our  assumption.  For 
more  complex  geometrical  and  topological  structures,  one  would,  however,  need  to 
compute  geodesic  distances  on  the  new  manifold  embedded  in  W. 

•  corrcoeff(*,  •):  is  the  linear  correlation  coefficient.  If  we  note  {dx}  and  {dz}  as  two 
ordered  sets  of  distances  (matrix  elements)  of  Dx  and  Dz,  respectively,  then: 

corrcoef(Dx,  T>z)  =  GXZ  ,  (8.3) 

axz  being  the  correlation  between  the  two  sets  {dx}  and  {dz},  and  ax  and  az  being 
the  standard  deviations  of  {dx}  and  {dz},  respectively. 

8.4.2  Simulation  examples 

We  saw  that  in  the  presence  of  noise,  the  performance  of  non-adaptive  isomap 
drastically  deteriorates  and  it  only  makes  sense  to  evaluate  performance  change  when  we 
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Figure  8.4:  Adaptive  embedding  of  the  noisy  swiss  roll  in  Figure  8.1. 


analyze  the  same  noisy  data  sets  using  the  adaptive  isomap.  We  consider  the  noisy  swiss 
roll  example  for  which  we  determine  an  embedding  as  shown  in  Figure  8.4.  We  next  treat 
a  slightly  more  challenging  case  constituted  of  two  adjacent  hemisphere  and  two  parallel 
sheets.  Figure  8.5  shows  the  results  obtained  for  the  hemispheres.  We  note  that  the  result  of 
the  non-adaptive  isomap  is  not  an  embedding.  This  is  due  to  the  connectivity  resulting  from 
using  a  Euclidean  neighborhood.  Indeed,  the  two  hemispheres  end  up  connected  through 
at  least  2  points.  We  avoid  this  connection  by  using  a  Mahalanobis  distance.  Figure  8.5  (c) 
shows  the  final  embedding  resulting  from  using  the  proposed  adaptive  isomap. 

In  Figure  8.6,  we  show  the  embedding  results  for  parallel  sheets  shown  in  Figure  8.6  (a). 


Figure  8.5:  Embedding  two  adjacent  hemispheres:  (a)  Hemispheres,  (b)  the  result  of  a 
non-adaptive  isomap  embedding,  (c)  the  result  of  the  proposed  adaptive  isomap. 


For  k  —  10,  we  see  in  Figure  8.6  that  non-adaptive  isomap  fails  again  to  define  an  embedding 
of  the  two  sheets.  The  reason  is  again  the  connection  that  occurs  when  using  a  Euclidean 
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distance  for  a  graph  construction.  The  result  of  the  adaptive  isomap  is  the  disconnection  of 
the  two  sheets.  As  they  should,  they  remain  two  distinct  structures  and  are  hence  separately 
embedded.  The  result  of  embedding  one  sheet  is  shown  in  Figure  8.6  (b).  Note  that  it  is 
exactly  the  sheet  itself,  but  now  in  2-dimensions  instead  of  3-dimensions.  For  a  different 
value  of  the  neighborhood  size  k  >  10,  we  find  the  results  in  (d)  and  (e)  for  the  adaptive 
and  adaptive  isomaps,  respectively.  We  notice  the  sensitivity  of  the  mapping  results  to  the 
neighborhood  size;  we  also  see  that  the  separation  between  the  two  sheets  did  not  happen 
in  both  cases;  however,  in  the  adaptive  case,  the  two  sheets  are  clearly  spread  on  a  plane, 
while  we  lose  one  sheet  in  the  non-adptive  case,  as  all  the  points  collapse  into  one  line. 


Figure  8.6:  Embedding  of  two  parallel  sheets:  (a)  the  initial  data,  (b)  Adaptive  isomap 
with  k  =  10  and  k  1  =  55.  Only  50%  of  the  initial  points  are  represented  here,  (c)  Non- 
adaptive  isomap  with  k  —  10.  100%  of  the  initial  points  are  represented  and  the  two  sheets 
are  overlapped.  Adaptive  isomap  (d)  and  non-adaptive  isomap  (e)  with  k  —  15.  100%  of 
the  initial  points  are  represented. 


To  further  evaluate  the  effect  of  noise  on  our  proposed  embedding  technique,  we 
increase  the  amount  of  Gaussian  noise  to  vary  between  0%  and  8%  of  the  orthogonal  distance 
between  the  two  parallel  sheets,  the  normal  distance  between  two  consecutive  levels  of  the 
swiss  roll,  and  the  orthogonal  distance  between  the  poles  of  the  two  adjacent  hemispheres. 
By  way  of  Monte  Carlo  simulations  on  the  data  in  hand,  we  obtain  the  results  shown 
in  Figure  8.7  (a)  and  (b).  We  establish  that  the  adaptive  isomap  technique  consistently 
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outperforms  the  non-adaptive  isomap  technique. 


(a)  (b) 


Figure  8.7:  Monte  Carlo  simulations. 


8.5  Network  failure  detection 

A  natural  source  of  data  which  share  similar  characteristics  with  much  of  the  data 
explored  above  is  that  of  sensors  in  a  given  network  topology.  We  propose  to  fully  exploit 
the  adaptive  isomap  algorithm  for  a  statistical  analysis  of  networks  in  the  presence  of  at¬ 
tacks  or  intrusions  with  the  ultimate  objective  of  detecting  them  and  potentially  localizing 
them.  Using  the  correlation  between  the  different  nodes  of  a  network  is  a  way  to  unfold  the 
interdependencies  between  them. 

Our  approach  starts  by  abstracting  networks  from  their  physical  layered  description;  we 
view  a  network  as  a  set  of  features/measurments.  The  same  feature  space  may  house  mul¬ 
tiple  networks  that  are  not  necessarily  the  same.  Our  hypothesis  is  that  each  network  is  a 
smooth  topologically  homogeneous  manifold  embedded  in  the  network  feature  space.  We 
consider  each  data  point  on  a  manifold  as  a  node  in  a  network.  In  addition,  an  attack  or 
any  intrusion  is  assumed  to  have  the  effect  of  destroying  the  initial  topological  structure  of 
a  manifold.  Our  intrusion  detection  task  is  thus  equivalent  to  detecting  and  determining 
any  topological  changes  on  the  nodes  manifold. 

Before  considering  the  dimensionality  reduction  of  a  network  data,  it  is  important  to  adap- 
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Figure  8.8:  Sets  connectivity  versus  neighborhood  size.  By  progressively  varying  the  size 
of  the  neighborhood,  we  identify  and  count  the  parts  that  get  merged.  We  talk  about  the 
persistence  of  these  parts  (Figure  8.9). 


Figure  8.9:  Detection  of  different  events,  (a)  Two  events  detected  when  using  a  Euclidean 
distance,  (b)  Three  events  detected  when  using  a  Mahalanobis  distance. 


tively  define  the  connectivity  between  the  different  nodes  of  a  network.  In  the  present 
experiment,  the  connectivity  is  dependent  on  the  correlation  between  the  measurements 
collected  at  each  node.  Using  a  mahalanobis  distance  becomes,  thus,  crucial  as  it  is  inde¬ 
pendent  of  any  scaling  we  may  apply  to  a  specific  set  of  measures.  We  illustrate  this  idea  in 
Figure  8.8,  where  we  observe  all  the  connectivity  events  associated  with  the  three  provided 
sets.  By  evolving  the  size  of  the  neighborhood,  we  connect  the  different  data  points  and 
find  the  number  of  the  final  resulting  structures  (connected  sets).  We  determine,  by  using  a 
Euclidean  distance  that  there  are  two  detected  events,  while  we  discover  three  events  when 
using  a  Mahalanobis  distance  (Figure  8.9). 
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Figure  8.10:  Distinction  between  attacked  nodes  in  red  and  normal  ones  in  blue. 

Given  the  correlated  and  noisy  nature  of  network  data,  we  adopt  the  adaptive 
isomap  algorithm.  We  use  the  DARPA  network  data  available  online  from  the  Knowledge 
Discovery  and  Data  (KDD)  Mining  Cup  1999  website.  We  use  all  41  parameters  to  define 
each  node,  i.e.,  N  =  41.  We  initially  only  use  data  that  is  subjected  to  a  single  attack.  We 
find  the  embedding  result  of  the  adaptive  isomap  in  Figure  8.10. 

The  first  comment  is  that  there  is  a  clear  separation  between  the  safe  and  attacked 
nodes;  nevertheless,  we  need  to  verify  that  this  visualization  makes  sense  and  that  we  are  in 
presence  of  an  accurate  embedding.  To  that  end,  we  look  at  the  embedding  dimensionality 
to  verify  that  it  is  indeed  equal  to  three.  We  find  in  Figure  8.11  (a)  the  residual  variance 
corresponding  to  this  experiment.  We  may  read  that  “3”  is  the  the  accurate  embedding 
dimensionality  when  the  network  is  subjected  to  one  attack;  however,  in  the  presence  of 
multiple  attacks,  the  embedding  dimensionality  increases  and  becomes  difficult  to  visualize. 

8.6  Conclusion 

We  proposed  a  new  adaptive  embedding  algorithm  that  first  learns  from  the  cor¬ 
relations  between  the  neighboring  points.  We  showed  that  combining  an  adaptive  distance 
with  an  existing  embedding  algorithm  leads  to  better  embedding  results  with  a  higher  ro¬ 
bustness  to  noise.  We  illustrated  this  technique  on  the  isomap  algorithm,  while  conceptually 
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Figure  8.11:  Embedding  dimensionality  of  networks  under  different  attacks,  (a)  Under 
multiple  attacks,  (b)  Under  a  single  attack. 


compatible  with  any  manifold  learning  technique  that  relies  on  a  connectivity  graph  for  the 
initial  data.  Our  initial  experiments  on  network  data  show  promising  results  that  we  plan 
to  further  investigate. 
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Chapter  9 

Future  research  directions 


In  this  thesis,  we  have  presented  different  techniques  to  compactly  represent  manifold¬ 
like  datasets,  with  an  emphasis  on  3-dimensional  objects.  All  our  research  work  was  primar¬ 
ily  motivated  by  applications  in  object  classification,  recognition  and  identification.  Many 
fundamental  questions  arose  in  the  course  of  this  work,  such  as  specifying  the  sampling  rate 
of  a  height  function  defined  on  a  complex  3D  shape,  remain  open  research  questions  and 
constitute  a  very  rich  avenue  for  future  research. 

Another  very  promising  and  new  theme  of  research  directly  impacts  networking  applica¬ 
tions,  as  quickly  illustrated  in  Chapter  8,  and  attempts  to  exploit  novel  structures  of  data 
which  may  approximately  lie  on  nonlinear  curved  spaces  or  manifolds.  For  concreteness  and 
conciseness,  we  concentrate  on  and  elaborate  on  these  two  directions  of  research,  the  first 
one  with  more  of  a  fundamental  interest  than  the  second. 

9.1  Fundamentals  of  shape  modeling 

In  order  to  achieve  an  efficient  modeling  of  2  and  3D  shapes,  a  good  understanding 
of  their  geometry  is  required.  To  that  end,  we  proceed  to  detail  the  crucial  elements  for  an 
accurate  shape  representation.  Specifically, 

1.  Optimal  reduction  of  triangulated  meshes. 

2.  Efficient  sampling  of  2-dimensional  surfaces  via  level  curves. 


3.  Surface  reconstruction. 
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9.1.1  Optimal  reduction  of  triangulated  meshes 

In  practice,  a  surface  S  is  arbitrarily  closely  approximated  by  a  triangular  mesh 
consisting  of  a  certain  number  of  nodes  (points).  As  discussed  in  Chapter  3,  each  surface 
S  is  characterized  by  a  parameter  5ft,  referred  to  as  the  characteristic  resolution  of  S.  This 
parameter  effectively  defines  the  minimal  number  of  points  required  to  ’’correctly”  represent 
a  shape  S. 

The  role  that  5ft  plays  in  reducing  the  computational  cost  of  representing  a  shape  with 
minimal  quality  loss  is  pivotal.  While  this  approach  has  demonstrated  a  representation  effi¬ 
ciency  improvement  of  up  to  an  order  of  magnitude,  the  illustration  remains  statistical  and 
empirical  and  hence  theoretically  unsubstantiated,  hence  raising  the  fundamental  question 
of  its  source  and  overall  viability. 

The  reduction  of  5ft  of  a  shape  as  shown,  was  equivalent  to  smoothing  the  shape, 
hence  affecting  it.  Analogously,  a  resolution  of  a  flat  norm  corresponds  to  a  certain  scale 
of  curvature,  which  in  turn  is  directly  related  to  the  characteristic  resolution  5ft  and  thus  to 
the  flat  norm  in  R3.  A  more  robust  and  systematic  modeling  of  3D  shapes  (for  instance), 
would  likey  incude  the  multi-scale  nature  of  the  flat  norm. 

9.1.2  Efficient  sampling 

We  have  seen  in  Eq.  (4.1)  how  a  surface  is  represented  by  a  set  of  disjoint  closed 
curves.  Sampling  a  surface  is  then  reduced  to  keeping  a  subset  of  these  curves.  The  obvious, 
and  crucial  question  is  what  to  keep  and  what  to  leave?  To  the  best  of  our  knowledge,  such 
a  sampling  rate  K  is  only  empirically  chosen  [1,  28].  Finding  the  effective  sampling  rate 
K  that  ensures  the  extraction  of,  and  only  of,  the  information  needed  to  represent  and 
reconstruct  a  given  shape,  would  thus  be  revolutionary  in  the  automatization  of  level  set 
based  algorithms. 

We  may  start  approaching  the  problem  by  first  considering  surfaces  (or  subsur¬ 
faces)  of  cardinality  one1.  Let  us  revert  to  the  probability  space  (V,^4,  Psr)  as  defined  in 
Section  4.2.3.  For  a  uniformly  sampled  space,  it  is  easy  to  note  that  the  instantaneous  value 
Pg^(i)  is  proportional  to  the  length  of  the  level  curve  at  the  value  i  of  the  GGF.  We  thus 
propose  to  measure  the  variability  A P  —  \Ps^(i)  —  P$i(j)\  induced  by  the  monotonous  and 

1  Cardinality  of  a  surface  is  the  number  of  closed  level  curves  for  a  given  value  of  the  considered  Morse 
function  on  that  surface. 
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continuous  evolution  of  the  level  curves.  In  fact,  it  is  this  variability  that  will  guide  an 
efficient  sampling  of  a  surface.  A P  is  nothing  but  the  change  in  the  length  of  C%  and  Cj 
up  to  a  constant.  In  summary,  the  idea  is  to  determine  A P  as  the  minimal  variation  to  be 
detected.  Then,  we  need  to  sample  Tfo  at  every  A P  variation. 

Many  questions  may  follow  from  the  proposed  preliminary  direction  in  investigat¬ 
ing  an  efficient  sampling  rule  for  surfaces.  One  may  consider  achieving  an  adaptive  sampling 
where  the  criterion  A P  changes  as  needed.  An  even  more  challenging  problem  is  to  consider 
complex  surfaces  where  the  cardinality  is  no  longer  equal  to  one. 

9.1.3  Surface  reconstruction 

Once  we  determine  how  detailed  the  representation  of  a  3D  object  needs  to  be,  it 
becomes  possible  to  ensure  a  good  surface  reconstruction.  Thus,  the  inverse  problem  to  the 
one  of  Section  9.1.2  may  be  posed,  i.e.,  we  are  given  few  level  curves  and  we  want  to  get 
the  initial  shape  back. 

Considering  the  observed  closed  sampling  curves  as  points  in  a  hyperspace,  we  may 
write  them  in  the  form  of  a  matrix  X,  where  each  column  corresponds  to  one  curve.  At  this 
level,  the  problem  is  reduced  to  interpolating  this  space  curve.  In  order  to  further  simplify 
the  problem,  we  propose  to  use  linear  algebra  and  compute  the  SVD2  decomposition  of  X. 

X  =  u(sVt)=UP.  (9.1) 

Our  problem  is  thus  reduced  to  finding  the  interpolated  version  of  P,  say  Q.  The  interpo¬ 
lated  version  of  X,  say  Y,  will  follow  as: 

Y  =  UZQ,  (9.2) 

where  Z  is  a  diagonal  filtering  matrix.  This  simple  technique  is  proposed  for  surfaces  of 
cardinality  one.  Again,  there  is  a  need  to  generalize  this  approach  to  generic  shapes.  It  may 
also  be  of  interest  to  combine  recognition  and  reconstruction  in  order  to  efficiently  restore 
occluded  objects. 

2 SVD:  Singular  Value  Decomposition 
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9.2  Behavioral  Modeling 

We  have  shown  throughout  this  thesis,  how  to  exploit  the  same  topological-geometric 
techniques  in  modeling  various  data  sets  (e.g.,  3D  objects,  sensor  network  data).  Modeling 
the  human  behavior  is  one  more  possible  setting,  that  may  benefit  from  our  proposed  set 
of  tools.  The  ultimate  goal  of  this  research  is  to  implement  an  automatic  system  capable  of 
understanding  complex  dynamic  scenes.  One  may  start  by  studying  elementary  activities 
(i.e.,  behaviors)  such  as  walking,  running,  and  jumping.  For  such  a  modeling  problem,  it  is 
of  interest  to  exploit  both  spatial  and  temporal  information  captured  in  video  frames.  We 
thus  may  view  each  frame  as  a  visual/spatial  description  of  the  studied  action  at  an  instant 
t.  The  description  in  question  is  the  precise  placement  of  the  subject  with  respect  to  his/her 
environment,  which  may  first  be  limited  to  the  floor.  We  propose  to  extract  this  description 


Figure  9.1:  Modeling  of  the  action  of  walking. 

through  the  contour  of  the  subject  and  its  relative  location.  We,  then,  may  recreate  the 
dynamic  scene  by  interpolating  the  surface  of  all  the  contours  parameterized  by  time.  By  so 
doing,  we  construct  the  one  dimensional  manifold  Ai  of  the  scene.  Our  modeling  consists 
therefore  in  embedding  M  in  R3  thanks  to  an  implementation  of  the  Whitney  embedding 
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theorem  (Section  5.3). 

We  investigate  the  discriminative  capabilities  of  the  proposed  method  by  looking 
at  three  scenes;  two  of  which  illustrate  the  action  of  jumping  (Figure  9.2)  and  the  third 
one  being  the  action  of  walking  (Figure  9.1).  The  preliminary  results  comfort  the  idea  of 


Figure  9.2:  Modeling  of  the  action  of  jumping  for  different  sequences. 
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a  signature  assigned  to  each  elementary  behavior.  There  are  nevertheless  many  questions 
that  one  may  needs  to  answer  in  order  to  fully  investigate  the  effectiveness  of  the  proposed 
technique.  Below  is  a  non  exhaustive  list  of  these  important  questions: 

1.  How  accurate  is  it  to  only  consider  the  exterior  contours? 

2.  How  much  does  the  change  in  topology  affect  the  final  model? 

3.  How  many  frames  are  needed  to  represent  each  movement? 

4.  How  to  deal  with  complicated  configurations  where  severe  occlusion  and  noise  are 
involved? 

5.  How  to  combine  many  actions  in  one  scene? 
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APPENDIX  A 


If  a  3D  shape  is  subjected  to  a  similarity  transform,  then  each  point  V*  in 


defined  by  the  coordinates  yXpVpZjj  is  transformed  as  follows 

V'}  =  aQY)  +  K,  (.3) 

where  a  is  a  scaling  factor,  Q  is  a  (3  x  3)  unitary  matrix  ,  i.e QQT  =  I3,  and  K  is  a 
(3  x  1)  translation  vector.  Considering  the  construction  in  (??)  and  the  transform  in  (.3), 
each  point  in  R3M  is  transformed  by: 

V7  =  aQ-V,  +  K,  (.4) 


,  and  K 


The  question  of  concern  is  to  find  the  transform  applied  to  a  modeling  curve  when  the  3D 
shape  it  originates  from  is  subjected  to  a  similarity  transform.  If  P  is  defined  by  (??),  we 
may  define  P7  as  follows 

P7  —  argmin  —  min  ||  PT\I/7  1 1 2  • 

Since  the  new  secants  Vi/7,  i  =  1  ,  •••  ,L,  result  from  the  same  equation  (.4)  as  V7,  we 
conclude  that  P7  is  directly  related  to  P  through  a  permutation  (or  a  rotation)  matrix.  In 
other  words:  P'  =  QPU,  where  UUT  =  I3;  thus,  we  may  write: 

W'  =  P/T  •  V', 

=  (QPU)T-  (aQV  +  K), 

=  a  (PU)T  •  Qr  •  QV  +  (QPU)T  •  K, 

=  aUr-PrV  +  L, 

=  a\JT  •  W  +  L.  (.5) 


From  (.5),  we  conclude  that  if  a  shape  in  3D  is  subjected  to  a  similarity  transform  (transla¬ 
tion,  rotation,  scaling,  or  a  combination  thereof),  then  its  modeling  curve  in  3D  is  subjected 
to  the  same  group  of  transforms. 
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APPENDIX  B 


Let  T (£),  N(t)  and  B(t)  be  the  tangent,  normal  and  binormal  vectors,  respectively, 
of  the  Frenet  frame  at  the  time  index  t.  Applying  Graham-Schmidt  normalization  to  this 
frame,  we  define  the  orthonormal  basis  (ei,  e2,  e3)  such  that: 


||e2||  =  e2  •  e2  =  1; 

(.6) 

e2  •  ei  =  0; 

(•7) 

e3  =  ei  x  e2. 

(.8) 

Our  task  is  to  define  the  two  turning  angles  ar{t)  and  asit).  We  recall  that  ar{t)  is  the 
signed  angle  formed  by  the  two  tangent  vectors  T(t)  and  T(t  +  1),  and  a  sit)  is  the  signed 
angle  formed  by  the  two  binormal  vectors  B  (t)  and  B  (t  +  1).  We  start  by  assuming  that 
T (t)  and  T (t  +  1)  construct  a  plane,  i.e.,  we  assume  that  ar{t)  is  not  null.  We  furthermore 
restrict  olt (t)  to  be  in  ] 0,  |].  It  follows  that  1  >  sin (ar(t))  >  0  and  1  >  cos(<ax(£))  >  0.  We 
may  write: 


ei  =  T  (f); 

e2  =  -(T(f  +  1)  —  aT(f)),  with  a,  b  ^  0, 

(.9) 


where  a  =  cos (ax(t)),  and  b  —  sin (ax(t)).  The  frame  at  (t  +  1)  needs  to  verify  the  following 
conditions: 


ei  —  T(t); 

e3  =  sign(N(t)  •  e3); 

,  _  T(t)  x  T(t  +  1) 

63  ““  ||T(t)  x  T(t  +  1)||’ 


(.10) 


e2  =  0  •  ei  +  cos  (olt)  •  63  +  sin(a)  •  e2;  (-11) 

or  e2  =  .  7  t  •  (N(t)  -  cos(a)e3).  (.12) 

sin(Q') 

Solving  these  equations  leads  to  the  following  result: 


©  =  sign(T(i  +  1)  •  N (£)  —  cos(<a)T (t  +  1)  •  e3). 


(.13) 
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Since  sign(T(£  +  1)  •  e^)  —  sign(T(t  +  1)  •  N (t)  —  cos(<a)T2  •  e%)  —  sign(T(t  +  1)  •  N (£)),  then: 

©  =  sign(T(£  +  1)  •  N (£))  •  arccos(T(£)  •  T(t  +  1)).  (.14) 


